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ABSTRACT 


This  rsssarch  used  Huygsns-Frssnst  wave  optics  computer  simulations  to 
investigate  the  effects  of  turtuience  strength  and  inner  sarie  on  the 
normaHzed  irradiance  variance  and  coherence  length  of  electromagnetic  waves 
pr(H>dg«fing  throunh  a  htftHilent  atmosphere.  These  mvestigattons  (feveloped 
several  fyMeUri,.  .or  udity  of  propagation  simidations  employtoig  a  numerical, 
split-step,  Huygens-r^/esrw*  method,  end  within  these  gukJelmes.  oonsidwed  five 
types  of  turbuienca  spectrum  inner  scale;  (1)  zero  imter  scale.  (2)  Gaussian  irmer 
scale,  <3)  Hiirs  arnf  (4)  FrehUch's  visoous^onvective  erihancement  inner  scales, 
and  (5)  turbuienoe  spectrum  truncation  from  the  discrete  grid  representation.  The 
simi^ation  results  showed  that  the  normalized  irradiance  variarx^  generally 
decreased  {-30%)  below  the  zero  inner  scale  vaNies  in  the  Rytov  regime  with 
increasing  inner  scale  size,  but  increased  monotonicaly  in  the  saturation  regime, 
and  agreed  within  2%  of  the  Rytov-T^arsW  predictions  at  tow  turbulence  strengths. 
The  E-flekt  coherence  length  in  a  spattaty  cortfined  beam,  with  either  spherical  or 
ptwie  wave  diyergertoe  arxf  zero  irmer  scale,  foiowed  the  Rytov-Tatarski-Fried 
predicitons  in  die  Rytov  regime,  but  deparfed  from  the  theory  in  the  sabiration 
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I.  INTRODUCTION 


Turbulwica  in  a  ^nrtiflad  fluid  cautaa  random  inhomoganaiflas  in 
tamparatura  and  tha  indax  of  rafracflon  that  acattar  and  diffract  a  wava 
propagating  throu(^  such  a  madkim.  Analytical  parturt>ation  tachroquas 
davalopad  ovar  tha  past  thraa  dacadaa  cannot  account  (Or  tha  variations 
obsarvad  axparfrnantally  in  tha  ampWuda  and  phasa  dMrttxjtiorM  of  a 
profMgatkig  wava  whan  turtulanoa  lavals  ara  high  arxUor  propagation  paths  are 
long  Ovar  tha  past  flffaan  yaars.  numartui  wava  optics  codas  baaad  on  tha 
Huygar^Frawtal  prindpia  hava  baan  davalopad  to  address  these  situations 
This  research  extended  these  codas  to  induda  inner  scale  in  tha  spactnjm  of 
rafTactiva  index  fluctuations  and  to  axamina  tfta  ooharanoa  length  and  tha 
affects  of  inner  scaia  in  the  Rytov  ragima  (low  turtMianoa  strengths  andfor  shod 
propagation  paths)  and  in  tha  saturation  ragima  (high  turtulanoa  strengths 
arHl/br  long  propagation  paths) 

As  a  frii  airtariWva  anaivsfal  and  *^»*»*^  davakioad  auktoHnas  for 
validity  of  computer  timulalions  employing  Huygans-f  rasnal  propagation  ovar 
muffipia  steps  (spM-stap  method)  These  guidsinas  mduda 
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•  For  an  NxN  grid  and  propagation  diatanoa  L.  the  grid  elemant  size  should 

be 

•  The  maxirnum  turbulerioe  strerigths  C.’ and  propagation  distances  L  for 
spherical  wave  propagation  with  an  NxN  grid  should  satisfy 

^  0.1  where  P«  •  0.497  C.*  k"*  L"'*  . 

•  Equivalentiy.  the  E*fleld  coherence  ler^jlh  r,  should  satisfy 

r,  2  2.5  Ax.  where  r,  represents  Fried’s  coherence  length  (for 
Sf^MKical  waves) 

•  Use  2  30  phase  screensfsieps  for  each  propagation 

•  Use  2  30  propagation  realizations  to  get  represerrtattive  statistics 

•  Phase  saeens  require  tow  spatiai  frequency  correction  to  gain  5% 
jccuracy  in  normalized  irradianoe  variance  and  as  ntuch  as  30%  in 
coherence  length 

•  Half  Width  at  half  maximum  of  the  irtmospheric  MTF  erto  an  iterative  fit  r, 
provide  the  most  stable  parametortzattons  of  the  E*field  coherence  length 

•  Telltale  signs  of  aliasing  incfude  a  fine-grained  kradtonoe  pattern,  a  boxed 
perimeter  to  the  irradiance  poOem.  and  peaking  of  the  energy  toward  the 
cantor  of  the  comput^wn  grid 

These  investigations  also  examined  four  chotoes  of  E-ftold  sotooe  function  and 
refined  the  methods  of  normaized  irradtonce  variance  and  E-ftold  coherence 
length  calarfation 

These  simutotions  tocorporatod  Ave  types  of  lurbulenoe  spectrum  inner 
scale  (1)  zero  inner  scale.  (2)  Gaussian  inner  scale,  p)  HVs  and  (4)  FrehliCh's 
viscous-convective  enhancement  inner  scales,  and  (5)  turtuienoe  spectrum 
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truncation  from  the  cMsaete  grid  representation.  For  the  more  physicaity 
realistic  viscouS'ConvecUve  enharwement  scale,  the  computer  simulations 
provided  the  foHovving  results; 

•  The  Hitt  «id  Frehllch  pman^erizations  performed  almost  identicaiiy. 
giving  less  than  3%  cNfferenoe  in  normattzed  irradlwioe  variance  over  the 
Rytov  and  saU^tton  regimes. 

•  In  the  Rytov  regime,  ttw  normalized  irradiance  variance  fOr  an 
approximately  spherically  diverging  beam  increased  (~30%)  and  then 
decreased  («30%)  compared  to  the  zero  Inner  scale  values  as  the  inner 
scale  size  ifKreased.  and  the  E-lleld  coherence  length  r^  decreased 
slightiy  (~5%)  compared  to  the  zero  inner  scale  coherence  length. 

•  In  the  saturation  regime,  tocreasmg  inner  scale  size  gave  morK)tonicaiiy 
increasmg  rMrmattzed  irradiance  variance  artol  increased  the  coherence 
length  (up  to  ^50%)  compared  to  ttMi  zero  inner  scale  case. 

The  effect  of  the  Gaussian  toner  scale  on  normalized  irradiance  variance  ar>d 

cohererxto  length  was  also  tovestigatod 

These  investigations  examined  the  behavior  of  the  E-fteld  ooherer>oe 

length  for  E>fields  (beams)  thiM  were  constrained  to  laleral  extent  to  the  grid 

size  but  whose  divergenoe  approxtourted  that  of  spherical  waves,  plane  waves, 

and  intermediate  beam  waves  The  con^Mitar  stoniattons  showed  that: 

•  For  the  Rytov  regime,  toe  spherical  and  plane  wave  cMes  gave 
coherence  len^^  witoto  5%  of  toe  Rytov-Talarolti>Fried  pradictions. 

•  In  the  saturation  regime,  toe  E-ltald  coherence  lengtos  for  toe  spherical 
wave  approximalion  deaeased  *25%  below  toe  theory,  whtte  the 
osherance  lertgtos  for  toe  plane  wave  approximatfon  varied  within 
•5%/^15%  of  toe  theory 
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•  The  E-field  coherence  lengths  for  the  beam  wave  cases  were  spaced 
between  the  spherical  and  plane  wave  values  for  the  Rytov  regime,  but 
increased  toward  the  spherical  wave  values  in  the  saturation  regime. 

The  spatially  confined  beams  used  in  these  simulations  showed  a  departure  in 

behavior  of  the  E-fieid  coherence  length  in  the  saturation  regime  from  the 

predictions  of  first  order  perturbation  theory. 

The  organization  of  this  dissertation  generally  follows  the  preceding 

summary  of  major  points.  Chapter  II  summarizes  the  theory  of  wave  optics  and 

the  impiert  antation  of  the  Huygens-Fresnei  propagation  in  a  computer 

simulation,  it  then  provides  an  overview  of  the  Rytov-Tatarski  theory  for  the 

effect  of  atmospheric  turbulence,  including  inner  scale,  on  the  normalized 

irradiance  variance  and  E-fieid  coherence  lengtti  as  applied  to  spherical  wave 

propagation.  Chapter  III  discusses  the  detailed  choices  of  physical  and 

simulation  parameters  and  develops  the  guidelines  to  ensure  validity  of  the 

simulation.  Chapter  IV  presents  the  computer  simulation  results,  first  for  the 

inner  scale  effects  on  normalized  irradiance  variance,  then  for  the  behavior  of 

the  E-fiekj  coherence  length  in  the  Rytov  and  saturation  regimes,  and  for  the 

effect  of  inner  scale  on  coherefioe  length. 
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II.  BACKGROUND 


A.  CHAPTER  OVERVIEW 

This  chapter  summarizes  first  the  theory  of  electromagnetic  wave 
propagation  and  the  Huygens-Fresnel  solution  to  the  scalar  wave  equation  and 
recasts  the  Huygens-Fresnel  solution  into  a  form  utilizing  fast  Fourier 
transforms  easily  implemented  in  a  computer  simulation.  It  then  summarizes 
the  Rytov-Tatarski  solution  to  the  scalar  wave  equation  which  used  first  order 
perturbation  theory  and  statistical  techniques  to  describe  the  spatial  variations 
of  *he  irradiance  of  an  electromagnetic  wave  propagated  through  a  turbulent 
medium.  The  Kolmogorov  spectrum  of  refractive  index  fluctuations  is 
introduced  as  well  as  five  types  of  inner  scale  which  modify  the  high  spatial 
frequency  portion  of  this  spectrum.  The  chapter  concludes  by  describing 
Fried's  method  for  parameterizing  the  coherence  length  of  the  E-field  for  a 
spherically-diverging  wave  propagating  through  a  turbulent  medium. 

a  WAVE  PROPAGATION 

Maxwell's  equations  describe  of  the  propagation  of  electromagnetic 
waves  forough  a  turbulent  medium.  Assuming  a  locally  homogeneous, 
isotropic,  and  linear  medium. 


V  -ci?  •  p. 

CD 

(2) 

(2) 

d.  ,i.  LStil, 

9t 

CD 

wfMf*  h«tt«d  (*}  quantttiM  mpnwarti  v*clort  Al  th«  frtqusnciM  and  ft«M 
strengths  of  toiterest,  the  medium  may  be  assumed  to  have  zero  iocei  free 
chi^  density  p  aryl  zero  (or  negiigibie)  ooryluctivily  o  Laser  beams  in  the 
atmosphere  satisfy  these  oorylitioris  ExparyJing  the  lefl-haryl  side  of  Eq  (1) 
aryl  cornbmmg  the  curl  of  Eq  (3)  and  the  dme  dertvalive  of  Eq  (4)  gives  a 
vector  wave  equation  for  the  E-fleid 

♦  V(#  •  V  In  €)  -  i»  f  ^  •  0.  (•> 

af* 

The  mkkfle  term  represents  depoiertzetion  of  the  E  ftetd  end  is  neglig«>ty  ameM 
for  propagation  through  a  turtuiwit  Mmosphere  (TetarM,  1961,  Lewrenoe  end 
Strohbehn.  1970)  Neglectir>g  the  middto  M»rm  simpMes  the  vector  wave 
equabon  to 

.  0.  (•) 

af» 
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!>)•  Fr»tn«(>Kirchofr  dittncOon  theory  provKiM  an  approidmM  tcily 
solution  to  this  vector  wove  equMon  Foioeting  HecM  (1987).  the  iphericel 
wi^  solution  m 

Eiu) .  •-«  ->.  m 

Where  E,  is  a  constant  and  k  ■  2s/X  Using  Eq  (7)  w(0)  lha  time  dependence 
factored  out.  the  Kirchoff  integral  theorem  becomes 

Which  must  be  evalueted  over  a  surface  S  endosing  the  field  poM  P.  Figure 
(1)  iWirtrtfirtes  the  relettonship  between  the  dtstarwes  ^  and  n  if  the  wavelenqth 
X  «  (.  n  .  Eq  (8)  reduces  to  the  Fresnef-Kiithofr  dWiracUon  formula 

•  T-'  //.^  '"*>  **  *** 

Where  K(6)  represents  the  obliquity  factor 

integrating  Eq  (9)  over  the  half  sphere  S  (shown  in  Fig.  (1))  as  the 
radius  approaches  infinity,  the  electric  field  amplitude  »  appreciable  only  over  a 
Ande  apedt^  area  A  of  the  hemisphere's  1UM  side  Then  Eq  (9)  reduces  to  the 
more  famitiar  Huygeris-Fresriei  expression 
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Plpupi  1  Frwntt  Wrchoff  dHtracOon  gtotntry  showing  rsiottonship  btwsn 
souros.  spsrturs  A.  surfisoo  8.  and  fItW  point  P. 

A  ^ 

Ths  oofnpm  aparturt  haidion  ooncaina  tha  aouroa  aphahcal  wava  factor, 
a]ip(itOi^  T?w  E-llaM  at  a  point  P  it  a  function  of  tha  E-fiald  ovar  a  finila 
apartura  A. 
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The  Huygene>Fresnei  formulation  .  Eq.  (10),  serves  as  foe  basis  for 
calculating  foe  propagation  of  a  wave  through  turbulenoe.  This  formula  can 
now  be  cast  foto  a  form  easily  implemented  in  a  propagation  computer 
simulation.  Applying  foe  Cartesian  coordinate  system  of  Fig.  (2),  foe  Huygens- 
Fresnel  prfodpie  becomes. 

CM  -tIL  --j - ■  'fW 


Aperture 


Flgurs  2  Propagation  coordfoate  system  showing  relations  between  foe 
aperture  variabfos  and  foe  field  variables. 
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^  reprtstnis  the  position  vector  of  the  observation  point  P  in  the  x-y  plane  and 
p  represents  the  position  vector  of  the  radiating  point  in  the  aperture  plane. 

The  paraxial  approximation  assumes  that 

\f\*z  ma  ip|<r.  (12) 

where  z  measures  the  lortgitiKflnai  distance  of  the  point  P  from  the  aperttfe. 
Consequently,  the  obliquity  factor  K<O)«(l4ooa0)/2  becomAS  approximately  1 
and  the  dervominator  of  the  integrarKt  is  approximately  z  .x>dman.1968), 

em  dA.  11») 


Fottowing  Roberts  (1986),  the  Fresnel  approximation  assumes  that  the  lateral 
displacement  between  the  radiating  point  and  the  observation  point 
I  p  -  p  I  «  z.  Then,  using  the  lowest  order  terms  for  the  series  expwision 
of  tfM  complex  expor)entii^s  argument,  Eq.  (13)  becomes 

.  -  -itefji  -  |iZf 

EVA  •  E{pSf)  •  ‘  *  ' 


XI 


fdpEipfli 


(14) 


.  a«y 

'  f;  ‘  I* 
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The  fBtilor  exp<-2xiz/X.)  represents  the  change  in  phase  of  the  optical  wave  from 
the  center  of  the  aperture  plane  to  the  plane  of  the  observation  point  P.  It 
applies  to  the  v^e  E-fieid  and  depends  only  on  the  propagation  distance  z.  It 
does  not  affect  the  phase  variation  across  the  E-field  nor  the  amplitude  of  the 
E-fiekj.  and  will  be  dropped  in  the  foltowirtg  expressions  for  simplicity. 

Equation  (14)  can  be  recast  into  a  form  utittzing  Fourier  transforms. 
Expandirni  the  aperture  field  E(^0)  in  Eq.  (14)  using  the  Fourier  transform 
identity 

£(p.O)  .  E(p',0),  0«) 

Eq.  (14)  becomes 

E(f,2)-  /df  /dp'  E(p',0)  ]  <“*•) 

(Note;  these  equations  use  spatial  frequency  /  (m  ')  instead  of  spatial 

wavenumber  k  >  2  s  /  (rad/m)  because  the  discrete  Fourier  transform 
computer  subroutines  used  in  these  investigations  are  written  with  spatial 
frequency  f  .)  Making  the  change  of  variables  f'mf-p  gives 
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E{t^  -  J(-dr)  [  /dip'  ]  • 


—1^1* 


(17) 


/df  •••''/dip'  •  E(p'.0)  [  fdt'  9 


-m*y  ^  XX 


i*  iVi  • 


]■ 


The  last  intagiation  over  it  the  Fourier  transform  of  a  Gaussian  function, 


/ 


df/  9'**^  a*' 


~i^l* 


Ur 


(18) 


Substituting  p for  p',  the  E-fiekf  of  Eq.  (17)  becomes 

E(Ar)  «  Jdf  a  '-Mfl*  /^p  £(^,0)  #  (1*) 

Equation  (19)  may  be  symbolicaNy  written  as 

E(Ar)  •  iF71  E71E(p,0)]I,  (2®) 

where  FT  and  IFT  represent  the  two-dtonensionai  Fourier  transform  and  inverse 
Fourier  transform,  respectively.  Equation  (20)  expresses  ttie  E-fiefd  at  a 
propagation  cHstanoe  z  in  terms  of  Fourier  transforms  of  the  E>fleld  at  z  »  0 
represented  in  a  Cartesian  coordinate  system.  This  form  of  the  Huygens- 
Fresnel  prindpie  is  useful  for  a  emulation  that  subdivides  the  propagation  path 
into  a  sequence  of  short  Fresnel  propagation  steps. 
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C.  EFFECTS  OF  TURBULENCE  ON  THE  PROPAGATED  WAVE 

1.  Imdhinot  VarltiiM  and  Innar  Scala 

The  next  a^  in  modeling  the  propagation  of  an  E-field  through  a 
turbulent  medium  requires  a  statistical  characterization  of  the  turbulence  and  its 
effect  on  the  spatial  statistics  of  the  propagating  wave.  To  begin  understanding 
the  effect  of  turbulence  upon  the  wave,  consider  a  spherical  wave  incident  upon 
a  mertium  rarxiomly  varyirig  index  of  refraction,  as  shown  in  Fig.  (3).  For 
smaH  variations  in  the  index  of  refraction  about  the  mean  value,  scattering 
occurs  predominantiy  in  the  forward  direction  (Tatarski,  1961).  The  regions  of 
index  of  refraction  fluctuation  accelerate/retard  portions  of  the  wavefront  over 
short  propagation  distances  and  cause  variations  in  phase  across  the  reference 
spherical  wavefront.  Subsequent  (tiffraction  and  interference  then  create 
variations  in  irradiance  across  the  wavefront.  Compared  to  propagation  through 
zero  turbulence  (Fig.  (4)),  a  sphericaRy  diverging  beam  now  exhibits  significant 
variations  in  phase  stkI  amplitude  (irradiance)  (Fig.  (5)). 

A  statistical  description  of  the  refractive  index  fluctuations  is 
needed  to  calculate  the  E4ield  variations.  The  index  of  refraction,  n,  in  air 
depends  upon  the  temperature  T  and  pressure  P  (Tatarski.  1961).  At  visible 
wavelengths, 
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Flgurt  3  Wav*  propagating  through  a  medium  having  random  iNK>mog6naitia8 
in  index  of  refraction. 


n  -  1  -  79x10^  (21) 

T{Ki 

The  temperature  profile  in  the  atmosphere  possesses  significant  stratification, 
as  shown  in  Fig.  (6)  (Waiters,  1994).  Velocity  shear  between  dHTerent  layers  in 
the  atmosphere  causes  turbulenoe  which  disrupts  the  interface  between  these 
layers,  mixing  regions  of  dMTerent  temperatures.  The  velocity  fluctuations 
generated  by  the  h^bulenoe  at  the  intertaoe  between  layers  generally  foNow  the 
Kolmogorov  spectrum. 
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Rgura  4  Irradiance  plot  for  a  spharfoa^  diverging  beam  propagated  through 
zero  turbulence. 


ngurof  Irradiano*  ptot  for  sph«rfcaly  diverging  bMm  propagate 
turtouiafKa. 


Ftgui*  •  Vertical  proflia  of  Itmparaliirt  in  the  atonotpharo. 

♦<«)  - 

aaauming  iaolropic  and  homoganaoua  turtuianoa  k  rapraaanta  tpattal 
wavanum{>ar  in  rad^  Tamparalura  iaa  Ipaaafva  additiva*  (Tatardd.  1961)  in 
ttia  atinoapfMca.  meaning  that  K  doaa  not  affaol  the  dynantica  of  the  turbulanoa. 
Thua.  the  thraa-dimanaional  tpadrum  of  tamparatura  fluctuationa  alao  fotioara 
the  Kolmogorov  apactnjm  Since  Eq.  (21)  raMaa  the  index  of  fafraction  to 
lamparaiura,  «ia  ipacarum  or  raaracava  moax  nucaiaDona  la  aiao  isovnogorov 
(TatariW.  1661) 

♦.(x,/)  -  atm  pfti)  x*”^.  (») 


17 


C,’(z)  is  a  cx>nsUint  of  profXM^onaNty  indicating  the  «drengtti  of  the 
tuftHiience.  in  reality,  this  spectrum  only  appUes  to  an  intermediate  range  of 
wavenumbers  known  as  dte  irterti^  subrange,  shown  in  Fig.  (7).  The 
boundaries  at  the  high  and  low  spatial  wavenumbers  are  the  inner  and  outer 
scales,  respectively,  and  are  cNscussed  below. 


Figure  7  Atmospheric  spectrum  of  retractive  index  fluctuations  showmg  inertiiri 
subrange  and  outer  and  Inner  scales 

\Mth  INS  urKterstandmg  of  Pie  effect  of  turbulence  on  a 
propagating  E-flald.  TatarsM  (1961)  derived  expressions  for  the  statistical 
proparHas  of  the  ampWude.  phase,  and  irradlanoe  of  a  spherical  wave 
propagating  through  a  turbulent  medkan  as  foiows  Consider  a  scalar 
component  of  the  vector  weve  equation  Eq  (6)  and  assume  a  hamtonic  time 
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dependence  for  the  E-field 

E{fJ)  -  u{n  (24) 

Perform  the  time  derivatives  and  introduce  the  wavenumber  k  =  2iUX  and  the 

irKtex  of  refraction 

/i  =  Cv^  (28) 

(where  c  »  velocity  of  light  in  vacuum,  p  »  magnetic  permeability  of  the 

medium,  e  -  electric  permittivity  of  the  medium  )  to  get  the  scalar  wave 
equation 

0.  (28) 

Expressing  the  index  of  refraction  as 

n{t)  *  1  f  riXt),  (27) 

where  ln,(A)  I  «  1,  and  usirtg  the  Rytov  metfrod  of  smooth  perturbations  that 
employs  u  *  exp(H0.  4'  *  V,  T,  +  ...,  and  u  =  u^  ♦  u,  +  ...,  the  wave 
equation  becomes 

♦  2ilr*n,(/)  -  0.  (28) 

This  has  a  sokrtion 

▼i(')  -  fhin  dV'.  (29) 
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for  a  tpharical  wav* 

u^\f\)  «  {O  .  oontfanfi  (30) 

In 


and  for  small  wavaiengtht  X  c  whara  ^  (tha  innar  scaia)  charactarizas  tha 
siza  of  tha  smaHast  fkictuatiorw  in  tha  tndax  of  rafraction.  tha  FrawMN 
approximation  in  Cartasian  ooordinatas  raduoas  Eq.  (29)  to 


aKp(  Mr 


zHx^*  ♦  (**♦/*)  -  2zz* 

_ _ 

z' (z-z^ 


Tha  statistics  of  tha  turtHjIanca  appaar  in  H'f  and  n,  in  tarms  of  spactral 
axpansions  whara  n,  oomairts  tha  thraa*dimansional  spatial  fraquancy  spactrum 
of  rafractiva  irtdax  fluctuations,  0,(x).  Tha  Ryfov  approximation  mtroduoas  tha 
log  ampWuda  fhjctuafion  X 

X  ■  m  ,  (S*) 

whara  A  and  A,  are  tha  ampiitudas  of  tha  turt)ular>ca  paftuft>ad  E-flald  u  and 
tha  firaa  spaca  E*fiald  u, 

"  ■  (Ml 

Uq  ■  Aq  • 

Than,  after  an  axtendad  sarias  of  man^ulations  (saa  Tatarski,  1961)  that 
aasuma  loctf  homogariaity  arxj  isotropy,  tha  varianca  of  tha  log  amplituda 
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fluctuatkm  X  of  a  spharical  wave  propagating  over  a  distance  L  becomes 


?  .  4****  //rtt  .  ♦.M)  <M) 

Where  k  •  2x/X.  and  k  «  spatial  wavenumber  (rad/m),  which  is  used  here 
instead  of  spatM  ftequency  f  (1/m)  to  be  consistent  with  previous  work. 
Assuming  that  the  irradiance  follows  a  log  normal  distribution  (Tatarski,  1961), 
the  normalized  irradiance  variance  is  a  function  of  the  log  amplitude  variance 

^  •  mp(4p )  -  1. 

The  assumption  of  isotropic.  homoger>eous  turbuierKe  gives  a 
Kolmogorov  spectrum  of  refractive  index  fluctuations  of  the  form 
^.(k.z)  *  0.033  C.’(z)  k'"'’  (Tatarski.  1961).  Assuming,  for  computational 
convenience,  that  this  k'"^  power  law  dependence  holds  over  the  entire  range 
of  spatial  wavenumbers  k.  and  that  the  turbulence  strength  C.’(z)  is  uniform 
along  the  path,  then  integration  of  Eq.  (34)  gives  the  log  amplitude  variance  for 
a  spherical  wave  as 

X*  -  0.124  Cj  (“» 

The  normalized  irradiance  variance  of  Eq.  (35)  becomes 
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of 


7* 


«(p(4  ?)  -  1  «  mpC  0^  Cj  -  1. 


(37) 


Low  turbutonce  (smaN  and/or  short  path  lengths  make  the  exponent  small 
enough  to  apply  the  approximation  exp(a) «  1*Kx.  giving 

^  -  0.407  Cf  -  pI  ,  (33) 

where  the  new  parameter  pj,  defirted  by  this  equation,  serves  as  the  baseline 
normalized  irradUmoe  variance  for  comparing  modifications  to  the  spectrum  of 
refractive  if>dex  fluctuations,  pj  also  ^iitates  plotting  normalized  irradiance 
variance  over  a  broad  range  of  integrated  turbulence  strengths. 

The  above  steps  characterize  the  effects  of  the  turbulent  medium 
upon  the  irradiar>oe  statistics  of  a  spherical  wave  assuming  a  simple 
Koknogorov  spectrum  of  refractive  ir>dex  fluctuatiorvi.  Incorporating  a  high 
spatial  frequency  roHoff  at  the  inr>er  scale,  as  shown  in  Fig.  (7).  modifies  the 
irradiance  statistics.  The  inner  scale  ^  represents  the  physical  size  where 
viscosity  of  the  medium  smooths  the  v^odty  fluctuations  by  dissipating  kinetic 
erHtrgy  and  tttermai  diffusion  smooths  the  temperature  fluctuations  (thus 
refractive  index  fiuctuattons),  removing  the  turbulent  character  of  the  medium  at 
this  scale.  This  smoothing  causes  a  rolioff  of  the  high  spatial  frequency  energy 
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of  the  refractive  index  spectnim.  Flatt6.  vytong.  and  Martin  (1993)  write  the 
three^imensional  isotropic  spectrum  as 


•(k,^)  »  0.033  F{kI^),  (5») 

where  F(k^  represents  a  particular  functional  form  of  the  inner  scale.  The 
parameter  consists  of  the  spatial  wavenumber  k  (radians  /  meter)  and  the 
inner  scale  parameter  (meters).  For  zero  inner  scale  (i.e.  no  high  spatial 
frequency  rollofr). 

=  1.  0  <  K  <  -.  (40) 

For  theoretical  and  computational  convenience,  a  Gaussian  rolloff  inner  scale  is 
often  used  (Tatarski,  1961)  to  represent  the  high  spatial  frequency  rolloff  from 
viscosity 


The  more  realistic  viscous-convective  enhancement  inner  scale  advocated  by 
Hill  (1978)  and  Frehlich  (1992)  exhibits  an  enhanced  spectrum  for  the 
waverujmbers  slightly  less  than  the  rolloff  point.  Viscosity  attenuates  the  kinetic 
energy  and  lowers  the  velocity  fluctuations  over  regions  near  the  inner  scale 
size  and  smaller  before  thermal  diffusion  can  smooth  all  the  temperature 
fluctuations  within  tiiese  regions.  This  disparity  alters  the  temperature  and 
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corrasponding  index  of  refraction  specfra.  Inner-8cale>8ized  patches  of  air  have 
different  temperatures  but  have  little  internal  velocity  variation.  Since  the 
Kolmogorov  spectrum  is  based  upon  the  spectrum  of  velocity  fluctuations,  these 
residual  temperature  fluctuations  cause  an  enhancement  to  the  Kolmogorov 
spectrum  around  the  inner  scale  size.  At  higher  spatial  frequencies  beyond  the 
inner  scale,  thermal  diffusion  does  smooth  out  these  residual  temperature 
fluctuations  and  the  spectrum  falls  off  sharply. 

Hill  provides  a  plot  of  F(k^  versus  for  the  viscous-convective 
enhancement,  that  is  suitable  for  implementing  the  viscous-convective 
enhancement  inner  scale  for  numerical  integrations  and  computer  simulations 
(Flattd,  Wang,  and  Martin,  1993).  Frehiich  presents  a  similar  characterization  of 
the  viscous-convective  enhancement  inner  scale  based  upon  laser  scintillation 
measurements  and  provides  a  four  parameter  fit  to  describe  this  version  of  the 
viscous-convective  enhancement  inner  scale. 

Since  most  numerical  simulations,  such  as  the  one  used  here, 
utilize  discrete  Fourier  transforms,  the  spatial  frequency  grid  mesh  chosen  for 
calculations  introduces  a  maximum  ^atial  frequency  established  by  the 
Nyquist  sampling  criterion  (discussed  in  Chapter  III).  This  limit  creates  a  grid 
cutoff  inner  scale  at  k,^. 


0  <  K  <  K, 
K  >  K, 


(42) 
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Figure  (8)  illustrates  these  five  inner  scales  for  eight  inner  scale 
sizes.  Solid  curves  represent  the  Hilt  version  of  the  viscous-convective  inner 
scale,  dashed  curves  represent  the  Frehlich  version  of  the  viscous-convective 
inner  scale,  and  dotted  curves  represent  Gaussian  inner  scales.  The  solid  line 
with  the  box  shape  represents  the  numerical  grid  cutoff  inner  scale  at  = 
318  rad/m  for  a  1024x1024  mesh.  The  inner  scaie  values  of  2,  3,  4,  5,  6,  7, 
10,  and  15  cm  refer  to  stratospheric  propagation  over  200  km,  assuming  an 
optical  wavelength  of  500  nm. 

The  outer  scale  l^,  corresponding  to  the  spatiai  wavenumber 
-  2ii/Lo,  represents  the  upper  size  limit  to  which  the  Kolmogorov  spectrum 
applies.  For  scale  sizes  larger  than  Lq  (or  spatial  wavenumbers  less  than  k,„), 
the  refractive  Index  fluctuations  level  off  to  a  finite  value  as  k  approaches  zero. 
Physically,  such  a  limitation  exists  as  the  spatial  wavenumbers  approach  zero 
because  the  refractive  index  fluctuations  cannot  become  arbitrarily  large,  or 
equivalently,  the  energy  represented  by  the  spectrum  must  remain  finite 
(Tatarski,  1961).  Estimates  of  the  outer  scale  for  the  stratosphere  lie  in  the 
range  of  tens  to  hundreds  of  meters.  However,  the  outer  scale  proves  to  be 
much  more  problematic  to  include  in  the  spectrum  of  refractive  index 
fluctuations  because  the  atmosphere  is  anisotropic  at  these  large  sizes 
(Tatarski,  1961).  The  Von  Kdrmdn  spectrum  (Tatarski,  1961), 
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Wavenumber,  K  (rod/m) 


FIgurt  8  Inner  scale  function  F(k^  for  8ie  grid  cutoff  (box).  Gaussian 
(dotted).  Hill  (solid),  and  Frehlich  (dashed)  inner  scales,  plotted  for  1024x1024 
grid.  L  s  200  km.  and  X  »  500  nm. 


•» 


0.033  Cj 

*  I*-”/*  * 


(43) 


which  has  been  used  to  Incorpoiate  an  outer  scale  for  some  calculations,  is 
based  upon  computational  convenience  more  than  physical  understanding. 
Additionally  the  computer  simulation  results  were  compared  against  the  Rytov- 
Tatarski  predictions,  which  were  derived  without  an  outer  scale.  Therefore, 
these  investigations  did  not  incorporate  an  outer  scale. 

2.  Atmoapherfc  MTF  and  Coheience  Length 

Propagation  through  a  turbulent  medium  not  only  affects  the 
irradiance  statistics,  but  also  reduces  the  spatial  coherence  of  the  wave  as 
characterized  by  the  transverse  coherence  length.  Fried  (1966,  p.  1372-1379) 
derived  a  long  exposure  modulation  transfer  function  (MTF)  for  a  spherical 
wave  propagating  through  a  turbulent  atmosphere  and  used  this  to 
parameterize  the  E-field  transverse  coherence  length,  r^.  Following  the  method 
and  notation  of  Fried  (which  used  spatial  frequency  f  in  (  m'^ )).  consider  the 
spatial  Fourier  transform  of  the  intensity  of  an  E-field  propagated  through  the 
turbulent  medium  and  imaged  by  a  thin  lens 

x(f)  •  Bf  Of  1/(2)  (44) 

where  u(2)  represents  the  E-fleld  in  foe  image  plane  and  B  represents  a 
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normalization  constant.  Fried  caiis  this  the  "MTF  of  the  image-forming  optical 
system**,  assuming  a  unit  impulse  (point)  illumination  and  anticipating  the  fact 

that  the  ensemble  average  of  t(/).  which  incorporates  the  effects  of 

atmospheric  turbulence,  turns  out  to  be  real.  Utilizing  the  Fourier  transform 
property  of  a  thin  lens, 

uiJi)  ^A[d9  Ui^) 

where  A  is  another  normalization  constant  and  U(^)  represents  the  E-fieid  in 
the  plane  of  the  thin  lens  aperture,  the  MTF  becomes 

t(/)  ^  A^Bf  d»  -  Xf^)  U{^).  (^) 

Expressing  the  E-field  U(  P)  as  the  product  of  a  zero  turbulence  propagation 
part,  W(^) ,  and  an  atmosphere-induced  perturbation,  V(^) , 

U(9)  -  W{9)  V^(P)  -  W{^ 

W( represents  the  uniformly  illuminated  aperture  function.  I{9)  represents 
fluctuations  in  the  log  amplitude  of  the  E-field  and  ^  represents  phase 
fluctuations. 

Taking  the  ensemble  average  <>  of  Eq.  (46)  over  many 
realizations  and  substituting  in  Eq.  (47)  gives  the  long  exposure  MTF 
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<^ih)iongmpo$um  =  V^(P)  >.  W 

Fried  denotes  the  expectation  of  the  fluctuations  as  the  "atmosphere's  MTP, 

MTF^f)  =  <  V{^  -  Xflf)  V{9)  \  (4«) 

and  Hufhagel  and  Stanley  (1964)  call  this  the  average  mutual  coherence  factor. 
Substituting  from  Eq.  (47). 

(  V*{9  -  XRf)  V{9)  )  =  \  (#0) 

Using  the  fact  that  ^  and  I  are  Gaussian  random  variables.  Fried  shows  that 
the  atmosphere's  MTF  reduces  to 

•  .  -4  (61) 
<  -  XRf)  )  =  •  *  ,  '  ' 


where  D(XR  |/|)  represents  the  wave  structure  function  and  is  related  to  the 
phase  structure  function  and  the  log  amplitude  structure  function  D,  by 

D(/)  »  D,(r)  ♦  D^{i),  when  r*  1 9-  XRf\.  (62) 


The  log  amplitude  structure  function  D,  and  phase  structure  function  are 
defined  as 


o,W  ■  ( |l(^  -  ), 

o«W  ■  ( l*W  -  ♦(»-x/V)i*  >. 


(SS) 
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•nd  cteicribe  th«  expactad  variation  in  tha  log  ampHtuda  and  phasa. 
ratpactivaly.  at  Uto  points  a  distanca  XRf  apart.  Danoting  tha  MTF  of  tha 
apartura  function  W( as 

t,(/)  ■  A*B  f  cf^  W{P),  (M) 


and  using  tha  tact  that  tha  atmospharic  MTF  of  Eq.  (51)  is  indapandant  of  P , 


tha  long  axposura  MTF  of  Eq.  (48)  now  bacomas 


ocxaifD 

$ 


(56) 


or. 


-4  o(W 

a  * 


(65) 


Tha  long  axposura  MTF  is  raiatad  to  tha  mutual  ooharanca  of  tha 
E-fiald.  Tha  mutual  ooharanca  function  (MCF)  dascribas  tha  autocorrelation 
batwaan  tha  E-fiald  at  two  points  and  is  defined  as  (Goodman,  1985) 

MCF  -  (  t;'(f,.f,)  >,  («T) 

whara  U  raprasants  tha  oomplax  E-fiald.  For  these  investigations  t^^t,  and 
explicit  time  notation  will  be  dropped.  Asstxning  homogeneity,  tha  MCF 
becomes  tha  spatial  autooorralation  of  tha  E-fiald 
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UCF{t)  -  /  df'  U*{F)  U{r  *  f). 


(M) 


Taking  tha  ansetnbie  avarage  of  Eq.  (46)  shows  that  tha  long  axposura  MTF 
comas  from  tha  avataga  mutual  coharenoa  furKtion  (autocorrelation)  of  the 
E-fiaid  in  tha  aperture  of  tha  Ians 

(t(/)>  ^  A*B  {f  -  X«f)  U(^  )  •  ( (  (/•(©)  U{^*XRf) ) ).  («•> 


Rewriting  Eq.  (56), 


<  W(t')  W{F*/)  ) 


m 


Equation  (60)  is  suitcbla  fOr  investigating  tha  wave  structure  function  D(r)  sir>oa 
computer  simulations  provide  tha  E*fialds  U  and  W  on  tha  right  hand  side  of 
Eq.  (60). 

Fried  (1966,  P.1380'1384)  dartvad  tha  wave  structure  furKtion  for 
a  spherical  wave  based  on  tha  thraa'dimansional  spatial  frequency  spectrum  of 
rafractiva  index  fluctuations,  ^.(k.z). 

0(1)  .  6«**»  dt  //  [1  -  ♦.M  .  dc.  (*1) 

where  z  represents  tha  position  along  tha  optical  path  length,  0  <  z  <  L  .  Fried 
analytically  solves  Eq.  (61)  for  a  spherical  wave  and  simple  Kolmogorov 
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furbultnot  that  hM  an  indttc  of  rsiricllon  itruclur*  ftjnclion 


0.(#yr)  -  («l 

and  a  thraa-cHmansioniri  apattai  tpadnim  of  rafractiva  vviax  fluctuationa 

#(k.^  -  0X»3  Ci(i)  k  (•») 

Suhatftut^  Eq.  (63)  into  Eq.  (61)  and  intagrattig  ovar  apatial  fraquancy  k  givaa 
D(r)  »  JL01  itr*  dz  Ci(z)  (^.  (64) 


Aaaummg  C.*(z)  s  oonatant  along  tha  opticai  path  givaa 

D(/)  .  1.0W  !(•  d  Lr^.  (•») 


Thia  aquation  raiataa  tha  wava  atnjctura  fijoctlon  0<r)  to  phyalcai  paramatara 
k«2ii/X,  C.^  and  L  for  tha  propagation  through  turtMianoa.  FoOowIng  Friad  arKi 
daflnk)g  tha  oonatant  r. 


4./ _ _ 

'iXWt  Ir*  Cf  i. 


(66) 


Eq.  (66)  tiaoofTiaa 

D{/)  .  tM  (-JT-  (6^ 

Tha  parama^  r,  charactartiaa  tha  coharanoa  lan(^  of  tha  E-I)ald  bacauM  tha 
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resolution  ailowod  by  the  turbulent  atmosphere  is  0^  •>  X  /  Tq  (Fried,  1966, 
p.  1380-1384). 

Substituting  Eq.  (67)  into  Eq.  (60)  gives 

( ( wn  uv'*D  >  >  (68) 

<  W{f*r)  >  ’ 

and  now  dire^  relates  a  measure  of  the  coherence  length,  r^,  to  the  E-fields 
produced  by  computer  simulation.  Similarly,  substitLting  the  integral  form  of  the 
wave  ^ructure  function,  Eq.  (61).  into  Eq.  (60)  and  numerically  integrating  for 
spectra  ^.(k,z)  with  different  inner  scales  allows  comparing  theory  and 
computer  simulation.  This  comparison  tadlitates  validation  of  the  computer 
simulations  incorporating  an  Inner  scale  at  low  turbulence  strengths  where  the 
perturbation-based  theory  remains  valid,  and  provides  a  mechanism  to  explore 
conditions  of  high  turbulence  strer>gth  and/or  k)r)g  propagation  path  lengths  (the 
saturation  regime)  where  the  perturbation  theory  may  no  longer  hold. 
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III.  COMPUTER  SIMULATION 


A.  PROPAGATION  CODE 

These  investigations  utilized  a  wave  optics  computer  simulation  code 
written  by  Brent  Ellerbroek  at  the  United  States  Air  Force  Phillips  Laboratory  in 
Albuquerque.  New  Mexico  that  incorporates  phase  screen  generation  routines 
written  by  Greg  Cochrane,  also  with  the  Riillips  Laboratory.  The  code,  known 
as  YAPS  (Yet  Another  Propagation  Simulation),  is  a  general  purpose  adaptive 
optics  simulation  code  written  in  FORTRAN  that  models  optical  propagation 
through  a  turbulent  atmosphere,  sensing  of  the  wavefront  with  Hartmann 
elements  (Hudgin,  1977}  and  a  CCO  array,  optimized  phase  correction 
calculation,  and  phase  compensation  via  a  deformable  mirror,  all  within  a  time 
indexed  framework.  These  investigations  utilized  only  the  propagation  portions 
of  the  code,  modified  the  code  to  parameterize  turbulence  strength  with 
instead  of  rg,  added  different  source  configurations,  incorporated  Gaussian,  Hill, 
and  Frehlich  inner  scales,  and  reduced  the  memory  requirements. 

The  computer  simulation  utilizes  the  sput-step  method  to  simulate 
propagation  of  the  E-field  through  tijrbulence,  as  illustrated  in  Fig.  (9).  A  source 
E -field  (left)  was  propagated  in  steps  (with  zero  turbulence)  over  the 
propagation  distance  L  with  a  random  phase  screen  applied  to  the  field  at  each 
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Z  -  0  Z  “  L 


Figuri  •  Computer  simulation  using  the  split-step  method.  The  source  (left)  is 
propagated  in  steps  out  to  z  »  L,  wWi  phase  screens  applied  at  each  step. 

step.  This  method  was  based  upon  the  extended  Huygens-Fresnei  principle 
(Yura.  1992) 

£,-jf Ea  -5^  Af(9)  •’  .  (W) 

which  is  Eq.  (10)  with  an  extra  e'^  factor  that  incorporates  the  random  variations 
in  log  amplitude  f  and  phase  ^ 

*»  .  •'  •'♦.  (70) 

For  paraxial  propagation,  if  the  distances  r\  are  short  enough,  diffraction  and 
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inMrfBTonos  not  havo  a  chanoa  to  cauaa  significant  log  ampiitijcle 
fiuctuations  f,  allowing  the  approximation 

e»  *  <7l) 

Inserting  this  approximation  into  Eqs.  (11)  •  (19)  gives 

E{r,z)  =  Jdp  E(^,0)  »/♦«  (72) 

Thus,  a  single  step  of  the  split-step  mefiiod  requires  applying  the  phase  screen 
e'*  to  the  E-field  and  propagating  the  result  a  distance  z  using  the  Huygens- 
Fresnel  propagator. 

To  implement  this  sequence  of  steps,  the  YAPS  code  followed  a  user- 
defined  list  of  tasks  (stored  as  an  input  file)  and  called  the  appropriate  routines 
to  accomplish  those  tasks.  For  these  investigations,  a  typical  list  (see  the 
Appendix)  first  set  up  the  initial  parameters,  including  random  number  seed, 
wavelength,  number  and  size  of  fields,  source  characteristics,  and  phase 
screen  size.  The  list  then  initialized  the  field  grid  and  applied  the  initial  source 
distribution,  propagated  the  field  in  steps,  generated  and  applied  phase  screens 
between  steps,  and  finally  saved  the  complex  values  of  the  propagated  field. 

The  YAPS  propagations  were  run  on  Sun  SparcStation-IO’s  having  128 
megabytes  of  RAM.  This  memory  consfi^int  allowed  a  maximum  grid  size  of 
1024x1024  (choosing  N  as  an  integer  power  of  2)  for  the  current  version  of 
YAPS.  Each  run  consisted  of  30  propagations,  each  with  32  phase 


screens/steps,  and  required  ~16  hours  to  complete.  These  investigations 
normally  utilized  nine  SparcStation-10's  simultaneously  doing  runs  with  different 
parameters.  In  all,  these  investigations  consumed  roughly  6000  hours  of 
computer  time.  A  small  Cray  mainframe  was  available  with  more  RAM,  but  was 
not  utilized  extensively  because  the  multiple  SparcStations  provided  more 
overall  computational  capacity. 

Once  the  YAPS  code  had  generated  and  saved  the  realizations  of  the 
E-field  propagated  through  turbulence,  separate  routines  written  in  interactive 
Data  Language  (IDL)  analyzed  and  displayed  the  fields.  Analysis  included 
displaying  two-  and  three-dimensional  plots  of  the  intensity  field,  calculating  the 
dependence  of  intensity  and  normalized  irradiance  variance  on  radial  distance, 
calculating  the  normalized  irradiance  variance  over  the  central  portion  of  the 
field,  calculating  the  atmospheric  MTF  and  corresponding  coherence  length, 
and  calculating  Strehl  and  intensity  ratios  (defined  in  Section  E  below). 

The  wave  optics  computer  simulation  required  many  choices  to  model 
the  stratospheric  propagation  scenario  and  to  ensure  validity  of  the  simulation. 
These  choices  included  wavelength,  propagation  distance,  inner  scale,  grid  size 
N,  the  physical  size  of  each  grid  element  Ax,  the  maximum  strength  of 
turbulence  the  source  function  of  tiie  E-field,  the  number  of  phase 
screens,  low  spatial  frequency  corrections,  number  of  realizations,  and  methods 
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of  calculating  average  statistical  values.  The  foliowing  sections  address  ttiese 
and  other  choices  and  develop  guidelines  for  vaiidity  of  the  simulations. 

a  PHYSICAL  PARAMETERS 

Wave  optics  computer  simulation  of  the  E-field  requires  selection  of  the 
parameters  that  describe  the  physical  properties  of  the  propagation.  Usually  a 
specific  propagation  scenario  has  been  selected  for  modeling,  giving  the 
desired  wavelength  X,  wavenumber  k=2ic4X,  propagation  distance  L.  and  the 
range  of  turbulence  strengths  C„^.  The  stratospheric  propagation  scenario 
chosen  here  used  X  =  500  nm,  L  =  200  km,  and  =  [IxlO"*’,  1x10'*l  m"*”. 
From  these  physical  parameters,  other  useful  scaling  parameters  occur,  such 
as  the  Fresnel  wavenumber  (Martin  and  Flatt6,  1988) 

•  (i)’"  =  ^ 


the  Rytov-Tatarski  normalized  irradiance  variance  (Flatt6,  Wang,  and  Martin, 
1993) 

^  =  a  Cj  .  (M) 


where  a-1.23  for  plane  waves  and  a-0.497  for  spherical  waves,  and  Fried's 
coherence  length  for  spherical  wave  propagation  (Fried.  1966,  p. 1380-1 384) 
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(78) 


^0  = 


6.88 


m 


2.91  //  dz  diz)  (!)•' 

=  ( - (/br  constant  d). 

M.09  k^d 


The  inner  scale  of  the  turbulence,  often  is  not  known  exactly,  but  values  can 
be  estimated.  For  stratospheric  propagation,  estimates  of  the  inner  scale  are 
around  1  - 15  cm  (Beland,  1993).  Once  the  inner  scale  is  known,  the  Hill  and 
Frehlich  versions  of  the  viscous-convective  enhancement  inner  scales  (Hill  and 
Clifford,  1978)  (Frehlich,  1992)  were  implemented  with  the  parameter  , 
where  k  represents  spatial  frequency.  The  outer  scale  Lg  for  stratospheric 
propagation  lies  in  the  range  of  tens  to  hundreds  of  meters  but,  again,  due  to 
the  difficulty  in  parameterization,  was  not  included  in  these  simulations. 


C.  COMPUTATION  GRID 

The  E-field  was  represented  as  an  NxN  array  of  complex  numbers  in  the 
plane  perpendicular  to  the  axis  of  propagation  .  Generally,  N  was  chosen  as 
large  as  possible,  consistent  with  the  amount  of  computer  memory  available 
and  the  time  required  to  run  a  simulation.  Larger  grid  size  N  allowed  a  wider 
spatial  extent  of  the  field  and/or  sampling  of  higher  spatial  frequencies  (denser 
mesh).  The  physical  size  that  each  grid  element  represented  had  to  be  chosen 
to  ensure  validity  of  the  computer  simulation. 
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Knepp  (1983)  discussed  the  relation  of  grid  element  size  Ax  and  physical 
grid  extent  NAx  to  inner  and  outer  scales,  to  proper  sampling  of  the  phase 
screens,  to  angular  spreading  of  the  fietd,  and  to  proper  sampling  of  the  spatial 
frequency  quadratic  phase  factor  in  the  Fourier  transform  formulation  of  the 
Huygens-Fresnel  propagation.  For  the  latter,  the  quadratic  phase  factor  in 
Eq.  (20)  takes  the  form 


♦ 


A  a.  * 


(78) 


where  k  is  a  spatial  wavenumber  in  rad/m.  Applying  the  Nyquist  criterion, 
which  requires  that  this  quadratic  phase  change  by  less  than  it  across  one  grid 
element  Ax,  Knepp  derives 

^  ^  2  NJ^xf  (77, 


Roberts  (1986)  applied  similar  techniques  and  derived  Eq.  (77)  without 
the  factor  of  2.  Again  applying  the  Nyquist  criterion  to  the  quadratic  phase 
factor  in  the  Fourier  transform  propagator  Eq.  (20)  but  utilizing  more  succinct 
differential  notation. 
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and, 


•  Af  2  It  X  L  4hw  < 


m 


where  represents  the  spatial  frequency  with  ttte  maximum  rate  of  change  of 
phase.  At  the  edge  of  the  grid,  f,^  is 


2  2  N  Ax  2  Ax 


(80) 


After  substituting  and  rearranging, 

L  <  (81) 


When  the  propagation  distance  L  is  known,  as  for  a  specific  propagation 
scenario,  and  when  the  grid  dimension,  N,  is  known,  then  the  grid  element  size 
is 

AX>^.  (M) 


Sampling  of  the  phase  factor  in  the  spatial  frequency  domain  must  meet 
the  Nyquist  criterion  as  a  minimum.  Some  circumstances  may  warrant  applying 
even  stricter  sampling  criteria  such  as  restricting  the  phase  to  change  by  less 
than  an  across  one  grid  element,  0  <  a  <  1.  Furthermore,  suppose  the 
propagation  involves  spatial  frequencies  out  to  where  0  <  p  <  1 .  Then 
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64  <  ax.  and  tha  maximum  apatiai  firaquancy  ia  Substituting  thasa 
constraints  into  tha  analysis  givas 


Ajr> 


n 

N 


(63) 


Whichavar  is  usad,  Eq.  (82)  and  (63)  provide  simple  formulas  for  choosing  the 
minimum  grid  alamant  size  for  tha  propagation. 

Spherical  wave  propagation  places  an  additional  constraint  on  the  choice 
of  grid  element  size.  In  the  parabolic  approximation,  a  spherical  wave  has  a 
quadratic  phase  curvature  which  represents  divergence  from  a  point  source  a 
distance  S  (focal  distance)  away 

4  *  4^' 

A  S 


where  p  measures  the  radial  distance  from  the  propagation  axis.  Roberts  (1986) 
analyzed  the  sampling  criteria  for  this  phase  just  as  for  the  spatial  frequency 
phase  (though  applied  to  the  case  of  a  two-step  Fourier  transform  Huygens- 
Fresnel  propagator  and  not  applied  specifically  to  spherical  waves).  Proceeding 
as  above  arfo  applying  the  Nyquist  criterion  to  this  phase  factor 


it .  A.  I  .  it. 

tfp  ftp  X  S  Ap 


(88) 
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(86) 


Ap  2  X  < 

A  5 


X. 


represents  the  maximum  radial  distance.  Using  corresponding 
to  the  edge  of  the  grid  and  Ap  ■  Ax. 


A 

PmtK  =  j  Ap  = 


Af  AX 
2  ‘ 


(87) 


Substituting  Eq.  (87)  into  Eq.  (86)  and  rearranging, 

4jr<|0.  (8») 


Again  generalize  the  analysis  by  requiring  the  maximum  phase  change 
across  one  element  to  be  ax  (0  <  a  <  1),  and  requiring  the  field  energy  to  be 
confined  within  a  region  of  radius  yp.^  ,  (0  <  y  <  1).  This  gives 


Ax< 


N 


XT7 

N  Y* 


(89) 


Combining  Eq.  (83)  and  (89), 

TTe 

N  Y* 


HI 

N  a 


<  Ax< 


(90) 


Spatial  frequency  sampling  considerations  for  the  Huygens-Fresnel  propagator 
have  placed  a  lower  limit  on  Ax,  while  spatial  sampling  considerations  of  the 
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quadratic  approximation  for  spharical  wava  phasa  hava  placad  an  uppar  limit  on 
Ax. 


Tha  computar  simulations  usad  in  thasa  invastigations  axaminad  high 
lavals  of  turbulanca  that  introducad  anargy  into  tha  highest  spatial  frequencies 
reprasantabla  on  the  grid  and  that  scattered  anargy  to  tha  edges  of  the  spatial 
grid.  Tha  simulations  also  assumed  that  tha  Nyquist  sampling  criterion  was 
sufficient.  These  considerations  spadfiad  tha  parameters  a  =  p  =  Y  ~  and 
lead  to 


<  Ax< 


(91) 


Additionally,  these  simulations  propagated  a  spherical  wave  from  the  point 
source  at  the  origin  (focus)  out  to  a  distance  S,  i.e.  S  =  L,  making  the 
inequalities  in  Eq.  (91)  become  the  equality 


(92) 


These  choices  unambiguously  determine  the  guideline  for  grid  element  size  for 
spherical  wave  propagation  based  upon  sampling  considerations  in  the  spatial 
and  spatial  frequency  domains.  Equation  (92)  also  determines  the  grid  element 
size  for  plane  wave  propagation  since  minimizing  the  grid  element  size  of 
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Eq.  (82)  optimizes  the  sampling  of  high  spatial  frequency  distortions  caused  by 
turbulence. 

Some  simplifications  were  made  in  the  above  analysis.  First,  the 
maximum  spatial  distance  and  maximum  spatial  frequency  used  in  the 
quadratic  phase  factors  were  chosen  for  the  nearest  edge  of  the  grid  and 
correspond  to  the  radius  of  the  largest  drde  that  will  fit  inside  the  grid.  These 
choices  disregarded  the  comers  of  the  grid,  but  this  omission  should  not  have 
affected  the  simulations  greatly  because  the  majority  of  the  energy  in  both  the 
spatial  and  spatial  frequency  domains  was  confined  within  the  radius  of  the 
drde  to  minimize  aliasing.  Second,  the  Nyquist  criterion  was  applied  to  the 
phase  change  across  the  x  or  y  dimension  of  the  element.  Analyzing  the  phase 

change  between  opposite  comers  of  a  grid  element  increases  Af  or  Ap  by  v/2  . 

Finally,  an  NxN  grid  with  N  even  requires  pladng  the  (x,y)  >  (0,0)  point  at  a  grid 
point,  such  as  (N/2  ■•'1,  N/2  ■•■I)  that  is  not  the  exact  geometric  center  of  the 
grid.  The  distance  to  the  nearest  side  is  (N-2)/2  elements  instead  of  N/2,  which 
is  a  minor  change  for  large  N.  Incorporating  tfiese  three  considerations  into  the 
analysis  for  a  spherical  wave  gives 
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How»v«r,  squaring  Eq.  (93)  wmi  raammging  laads  to 


i<i_«Us 


♦  (/V-2)‘ 

so  that  a  spherical  wave  propagation  with  these  additionai  constraints  should 


m 


only  be  carried  out  over  a  maximum  distance  less  ttian  S/4.  Since  these 
investigations  needed  to  propagate  a  spherical  wave  over  the  full  focus 
distance  S  starting  at  the  source,  the  simpler  expression  Eq.  (92)  was  used  to 
determine  the  grid  element  size  fOr  these  investigations.  As  a  comparison, 
FlattA,  Wang,  arni  Martin  (1903)  used  a  grid  eiement  size  of  0.7mm  fOr  a 
1024x1024  grid  and  0.5mm  for  a  2048x2048  grid  with  XL  =  0.000638  m^. 


Equation  (92)  with  these  parameters  prescribes  grid  element  sizes  of  0.78  mm 
and  0.55  mm,  respectively,  which  are  ~10%  larger  to  meet  sampling 
considerations  for  the  Huygens-Fresnei  propagator. 

The  spW-step  method  in  the  ccxnputer  simulation  divides  the  optical  path 
into  multiple  steps.  However,  the  disti^Yce  L  used  to  determine  the  grid 
eiement  size  must  correspond  to  the  total  propagation  path  length  and  not  the 
step  size.  To  justify  tNs,  consider  the  vacuum  propagation  of  a  field  across  die 
distance  L  »  Az.  If  prc^xagation  occurs  in  a  sin(^  step,  then  the  Fourier 
transform  Huygens-Fresnei  propagator,  Eq.  (20).  becomes 


E(Az)  -  /FTl  e  '**^'*  F71E(P)]]. 
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If  ttw  path  has  two  steps,  then 


e(Ar)  -  fn£<Az/2)ll. 


(H) 


and 


£(Az/2)  =  iFT[  9  »  FT\EiO)]]. 


(97) 


Substituting  Eq.  (97)  into  Eq.  (96), 


_j_«  1-1  Az«i  fQa\ 

E(Az)-/F71a  *  Frf/Fn*  *  PT[E(0)]\]\.  '  * 


But  ihe  middle  Fourier  transform/  inverse  Fourier  transform  pair  cancel  each 
other,  ar>d  this  two  step  vacuum  propagation  becomes  the  one  step  propagation 
of  Eq.  (95).  Thus,  the  grid  element  size  must  be  chosen  to  correspond  to  the 
total  distance  Az  -  L  since  multi-step  vacuum  propagations  mathematically 
reduce  to  a  single  step  of  L. 

D.  SOURCES 

The  choice  of  a  given  propagation  scenario  significantly  affects  the 
irradiance  and  coherence  statistics.  Plane  wave  propagation  differs  from 
spherical  wave  propagation,  as  seen  in  the  Rytov-Tatarski  theory  (Tatarski, 
1961),  while  beam  wave  propagation  (Gaussian  intensity  profile)  exhibits  an 
intermediate  behavior  (ishimaru,  1978).  17ie  source  must  also  be  chosen  to 
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keep  the  lateral  extent  of  the  source  propagation  within  tfie  spatial  and  spatial 
frequency  limits  of  the  computer  simulation. 

Coherent  plane,  beam,  and  spherical  waves  differ  in  the  shape  of  the 
isophase  surfaces  of  the  E-field.  Beam  wave  sources  possess  a  quadratic 
phase  (tshimaru,  1978), 

Mp)  ■  p*.  <»•) 

where  Rg  represents  the  radius  of  curvature  of  the  wavefront.  Plane  wave 
sources  have  an  infinite  radius  of  curvature  so  that 

^p)  *  constant.  C**®®) 

Spherical  wave  sources  have  a  radius  of  curvature  equal  to  the  propagation 
distance  from  the  origin,  or  focus.  S 

♦(p)  =  («’) 

Turbulence  introduces  random  phase  shifts  across  the  wavefront  while 
diffraction  and  interference  further  distort  the  wavefront  during  propagation 
causing  the  light  to  become  partially  coherent.  The  resulting  partially  coherent 
E-field  no  longer  possesses  a  simple  isophase  surface  and  the  coherent 
wavefront  characterizations  no  longer  apply.  The  differences  between  plane, 
beam,  and  spherical  propagation  appear  in  differences  in  statistical  properties. 
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such  as  normalized  irradiance  variance  (Tatarski,  1961) 


^2  ^J23C^  pkne 

\).487  Cj  spherical 


(102) 


(see  Ishimaru,  1978  for  the  more  complicated  beam  wave  expression),  and  the 
coherence  length  (Fried,  1966,  p.  1380-1384) 

3.02  (Cj)"^  k-^  spherical 
ro  =  {  (^03) 

1.68  k^  Z.-^,  plane. 

When  the  E-fieid  has  propagated  through  turbulence  to  the  far  field  of 
the  source,  diffraction  has  caused  the  E-field  to  expand  laterally  so  that  its 
statistical  properties  approach  those  of  a  spherical  wave.  The  far  field  begins 
at  a  distance  (Saleh  and  Teich,  1991) 

Z  -  10  (104) 


Where  p  is  the  radius  of  the  source.  Spherical  propagation  properties  result 
when  the  maximum  radius  of  the  source  approximately  satisfies 

\  z.  C”*®*) 

For  a  stratospheric  propagation  scenario  with  L~200  km  and  \  =  500  nm,  the 
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source  radius  must  be  less  than  10  cm  to  obtain  enough  divergence  to  prochjce 
spherical  propagation  statistics  at  z  ~  L. 

Plane  wave  propagation  statistical  properties  result  when  the  E-field  is 
evaluated  in  the  near  field  of  the  source  because  the  E-field  does  not  have  the 
opportunity  to  diffract  or  expand  significantly.  Starting  with  Eq.  (105)  for  the  far 
field  point,  near  field  propagation  satisfies 

z.±sL«  10^,  (10S) 

10  X  X 


or. 

>  ^id  X  z.  007) 

For  the  stratospheric  propagation  scenario,  source  radius  >  1  m  for  near 
field  propagation.  However,  validity  of  the  Fresnel  approximation  to  the 
Huygens-Fresnel  equation  places  an  upper  bound  on  p^  (Saleh  and  Teich, 
1991) 

(p5«)*  «  4  X,  (“108) 

(derived  from  considering  the  Taylor  series  expansion  of  the  Cartesian 
expression  for  r).  For  the  L  =  200  km,  X  =  500  nm  stratospheric  propagation 
scenario,  the  upper  bound  on  source  lateral  extent  becomes  <  350  m.  A 
1024x1024  grid  with  the  grid  element  size  Ax  =  0.99  cm  from  Eq.  (92)  gives  a 
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maximum  grid  radius  of  only  5  m,  so  tiiat  the  Fresnel  approximation  is  satisfied 
for  any  source  represented  on  this  grid. 

Ishimaru  (1978)  evaluated  the  irradiance  statistics  for  the  intermediate 
Gaussian  beam  wave  case.  The  normalized  irradiance  variance  along  the 
beam  axis  depends  upon  the  beam  waist  size,  Wg,  and  the  radius  of  curvature, 
Rq.  Numerical  integration  of  the  log  amplitude  variance  formula  (Walters,  1994) 
provides  a  smooth  transition  from  spherical  wave  variance  statistics  to  plane 
wave  statistics,  with  a  dip  in  between,  as  shown  in  Fig.  (10) .  Numerical 
simulation  results  with  Gaussian  and  Airy-type  sources  of  varying  widths  have  a 
corresponding  behavior,  as  shown  in  Fig.  (11). 

The  finite  grid  in  a  computer  simulation  places  limitations  on  the  source 
E-field.  The  physical  grid  element  size  sets  a  lower  bound  on  the  width  of  a 
narrow  source  approximating  a  point  source.  A  source  of  width  dose  to  a 
single  grid  element  may  still  be  undersampled,  causing  the  resulting  propagated 
field  to  exhibit  sidelobes  from  absence  of  the  proper  high  spatial  frequences  in 
the  representation.  Correspondingly,  a  spherical  wave  from  a  very  narrow 
source  will  propagate  with  large  angular  divergence  that  may  exceed  the 
physical  dimensions  of  the  grid  and  cause  energy  to  leak  off  the  grid  and  be 
aliased  back  into  the  field.  Thus,  the  source  must  be  wide  enough  to  constrain 
the  propagated  field  so  that  most  (>90%)  of  the  energy  remains  within  the  grid. 
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Normalizad  irrodionc*  vorionee  /  Bcto  squorad  (tpharieai) 


Q 


Figura  10  Normalized  irradiance  variance  versus  beam  waist  size  Wg,  wl 
Q^icWg’/XL,  from  numerical  integration  of  Ishimaru's  predictions  for  beam 
waves. 
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Normalized  irrodionce  vorionce  /  Beta  squared  (spherical) 


Figura  11  Normalized  irradiance  variance  versus  beam  waist  size  Wq,  where 
rbffWo^/XL.  Solid  line  =  Airy-type  source;  dashed  line  =  Gaussian  source. 
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Similwty.  a  plane  wave  source  cwmot  extend  too  near  the  grid  edge  because 
turbulence  can  scatter  energy  off  the  grid  only  to  be  aliased  back  in. 

Statistical  calculations  on  the  final  beam  pattern  require  a  reasonably 
uniform  central  patch  to  perform  compt«tation8.  If  the  region  over  which 
statistical  calculations  are  made  has  variations,  such  as  sidelobes  arising  simply 
from  vacuum  propagation,  these  variations  become  a  part  of  the  statistics,  such 
as  the  normalized  irradianoe  variance.  Often,  Gaussian  type  profiles  present  a 
minimum  of  such  variations  because  they  do  not  have  as  much  energy  in  high 
spatial  frequencies.  However,  because  they  are  rK>t  fiat  over  the  calculation 
region,  different  regions  of  the  beam  must  often  be  weighted  with  respect  to 
their  mean  intensity  in  doing  statisticai  calculations  such  as  the  normalized 
irradiance  variance. 

To  meet  these  constraints,  Martin  and  Flatt4  (1990)  applied  a  quadratic 
phase  curvature  to  a  Gaussian  source,  thereby  increasing  the  divergence  and 
flattening  the  final  irradiance  pattern.  Spherical  wave  sources  were  simply 
chosen  narrower  than  the  plane  wave  sources  to  make  them  diverge  more. 
Flatt6,  Wang,  and  Martin  (1993)  also  used  a  super-Gaussian  to  model  an 
extended  beam  source. 

The  simulations  presented  here  used  an  alternate  method  suggested  by 
Ellerbroek  (1993).  The  final  field  at  z  -  L  was  specified  as  an  aperture  of 
radius  equal  to  one  half  the  grid  radius  and  given  a  quadratic  phase 
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corresponding  to  a  spherical  wave  originating  at  the  origin  z  =  0  (Fig.  (12)). 

This  aperture  field  was  then  backward  propagated  without  turbulence  from  z  =  L 
to  the  focus  z  >  0,  effectively  Fourier  transforming  the  field  and  yielding  Vne 
numerical  equivalent  of  the  Airy  pattern  (Fig.  (13)).  This  source  was  then  used 
as  the  approximation  to  a  point  source  for  all  spherical  wave  propagations.  The 
advantages  of  this  method  were  that  the  source  was  quite  narrow,  having 
appreciable  amplitude  over  only  a  few  grid  elements,  and  that  the  zero 
turbulence  propagated  field  at  z  -  L  was  necessarily  constrained  within  the  grid 
and  possessed  a  central  region  with  uniform  illumination. 

The  significant  energy  at  high  spafiai  frequencies  required  to  represent 
the  sharp  edges  of  the  initial  aperture  at  z  -  L  created  one  drawback  since  it 
caused  strong  Fresnel  fluctuations  at  intermediate  distances  between  0  and  L. 
To  investigate  the  significance  of  these  Fresnel  fluctuations,  the  energy  at  these 
high  spatial  frequencies  were  reduced  by  windowing  the  aperture  before  the 
backward  propagation  by  applying  a  Gaussian  rolloff  to  the  cylinder  edges. 


1^ 


r<  0.7  r' 
r  z  0.7  r',  n=2, 


(109) 


where  r*  represented  the  corresponding  aperture  radius.  Because  this  rolloff 
reduced  the  energy  in  high  spatial  frequencies  for  the  source  representation, 
this  slightly  increased  the  size  of  the  final  central  region  of  the  E-fietd  where 
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Figure  12  Cross  section  of  intensity  profile  for  final  propagated  beam  showing 
confinement  of  energy  and  uniform  central  illumination  (beam  radius  =  70  grid 
elements  here). 


Figure  13  Cross  section  of  intensity  profile  of  numerical  equivalent  to  Airy 
source. 
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statistical  calculations  could  be  performed.  Other  values  for  the  power  n  on  the 
exponent  were  tried  (n  s  1,  n  s  5,  n  =  10)  but  the  simple  Gaussian  (n  =  2) 
worked  best.  While  this  rolloff  appeared  useful,  identical  runs,  one  with  the 
sharply  defined  aperture  edge  and  one  with  a  rolled-off  Gaussian  edge,  showed 
less  than  0.5%  difference  in  the  normalized  irradiance  variance  values.  This 
difference  was  negligible  compared  to  statistical  fluctuations  in  normalized 
irradiance  variance  of  up  to  ^10%  between  runs  ^th  different  random  phase 
screens.  Empirical  observation  of  simulations  revealed  that  once  even  modest 
amounts  of  turbulence  existed  (typically  ^  within  the  Rytov 

regime),  the  high  spatial  frequency  energy  introduced  by  the  turbulence 
dominated  the  details  of  the  source. 

Another  deficiency  was  the  representation  of  the  initial  cylindrical 
aperture  on  a  rectangular  grid.  Due  to  the  Cartesian  nature  of  the  grid,  the 
curved  beam  edge  was  actually  jagged  instead  of  circular.  However,  for  large 
enough  grid  (e.g.  1024x1024)  this  departure  from  a  cylinder  introduced 
negligible  effect.  The  only  case  where  a  difference  was  noted  was  with  a 
256x256  grid,  and  increasing  the  beam  radius  a  small  fraction  of  a  grid  element 
size  eliminated  that  difference. 

Another  limitation  of  this  method  was  that  the  final  E-field  could  not 
approach  the  grid  radius.  If  this  happened,  then  high  turbulence  strengths 
would  scatter  energy  off  the  grid,  producing  aliasing.  Due  to  the  periodic 
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Fourier  transform  method  of  propagation,  this  energy  would  not  be  lost,  but 
aliased  back  in  at  lower  spatial  frequencies.  Martin  and  Fiattd  (1990) 
implemented  an  attenuating  region  just  inside  the  grid  radius  to  absorb  this 
energy  and  to  prevent  its  being  aliased  back  in.  This  obviously  introduced  a 
type  of  windowing  function  in  the  spatial  domain  and  did  not  conserve  energy  in 
the  field,  possibly  complicating  final  field  comparisons.  Rather  than  include 
those  effects,  these  simulations  simply  chose  the  aperture  radius  suffidentiy 
small  (at  one  half  the  grid  radius)  to  prevent  significant  energy  scattering  off  the 
grid  in  the  spatial  domain  while  still  providing  a  relatively  uniformly  illuminated 
central  region  for  statistical  calculations. 

Obviously,  other  choices  for  a  source  existed.  For  example,  the  initial 
field  could  have  been  expressed  according  to  an  analytic  expression  for  the 
Airy  pattern  (Fig.  (14)).  This  source  gave  a  vacuum  propagated  final  field  that 
was  very  dose  to  a  broad  cylinder,  but  which  exhibited  noticeable  sidelobes  at 
the  edges  (Fig.  (15)).  These  sidelobes  came  from  spatial  truncation  of  the  Airy 
pattern.  The  back  propagated  numerical  Airy  pattern  of  Fig.  (13)  differs  from  an 
analytic  Airy  pattern  in  a  way  that  eliminates  the  sidelobes  in  the  final  irradiance 
field. 

A  point  source  could  be  modeled  by  a  single,  nonzero  pixel  at  the  center 
of  the  grid  being  nonzero,  as  shown  in  Fig.  (16).  Figure  (17)  shows  that  the 
vacuum  propagated  field  exhibited  noticeable  sidelobes  at  the  edges  and 
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covered  the  entire  aperture,  meaning  that  turbulence  would  immediately  scatter 
energy  out  of  the  grid  only  to  be  aliased  back  in  and  distort  the  field.  A  similar 
source  was  a  uniformly  illuminated  grid  at  z  >  L  with  quadratic  phase,  backward 
propagated  to  z  =  0.  While  this  eliminated  the  vacuum  propagation  ringing 
problem  by  definition,  it  immediately  suffered  from  aliasing  when  turbulence  was 
nonzero  and  scattered  energy  off  the  grid. 

E.  MAXIMUM  TURBULENCE  STRENGTH 

One  of  the  more  difficult  choices  in  a  computer  simulation  of  propagation 
is  determining  the  range  of  turbulent  strengths  over  which  a  simulation  is 
valid,  in  general,  the  processes  of  optical  propagation  through  a  turbulent 
medium  are  not  band  limited  in  spatial  frequency.  But  when  a  computer 
simulation  uses  a  finite  grid  to  implement  a  source  function,  to  propagate  via 
the  Huygens-Fresnei  prindpie,  and  to  model  the  turbulence  by  phase  screens, 
aliasing  inevitably  occurs  due  to  the  finite  sampling  interval.  This  problem  only 
becomes  more  severe  as  turbulence  strength  increases.  The  coherence  length 
measures  the  physical  distance  over  which  the  mutual  coherence  of  the  E-field, 

( E*  E  ),  dedines  to  e  ’  of  its  peak  value.  As  turbulence  increases,  the  E-field 
fluctuates  significantly  over  smaller  and  smaller  distances  and  the  coherence 
length  decreases.  An  NxN  discrete  grid  representation  of  the  E-field  samples  at 
a  specific  minimum  distance,  and  hence  at  a  maximum  spatial  frequency.  If  the 
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E-field  ^uctuates  in  less  than  this  minimum  distance,  the  grid  samples  will  not 
represent  the  E-field  accurately,  causing  aliasing.  The  question  is  not  whether 
aliasing  occurs,  but  how  much  aliasing  occurs  for  a  given  set  of  simulation 
parameters,  and  at  what  point  does  aliasing  invalidate  the  results  of  the 
simulation. 

These  investigations  identified  three  telltale  signs  of  aliasing.  First,  as 
aliasing  became  significant,  irradiance  plots  of  the  turbulent  E-field  lost  structure 
and  exhibited  an  isotropic  fine-grained  appearance,  as  Fig.  (18)  shows. 

Second,  the  irradiance  that  should  have  been  roughly  radially  symmetric 


Figure  18  Intensity  plot  showing  fine-grained  pattern  due  to  significant  aliasing. 
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became  bounded  within  a  fuzzy  rectangular  region,  as  Fig.  (19)  illustrates.  The 
irradiance  at  the  center  of  the  image  increased,  or  started  peaking,  while  the 
irradiance  at  the  edges  decreased,  as  shown  in  Fig.  (20). 

Figure  (21)  also  illustrates  this  peaking  behavior  by  plotting  average 
irradiance  versus  radius,  in  units  of  grid  element  size.  The  E-field  propagated 
with  zero  turbulence  had  the  radially  symmetric  aperture  profile  in  this 
simulation.  As  the  turbulence  became  stronger,  the  sharp  edges  of  the  initial 
cylindrical  irradiance  pattern  became  rounded  and  energy  spread  outward. 
However,  when  significant  aliasing  started  occurring  around  ~  5  for  this 
64x64  ghd  representation,  the  energy  began  creeping  inward  and  the  center  of 
the  field  started  increasing  in  irradiance. 

The  irradiance  at  the  edge  of  the  grid  was  -1/40  that  of  the  center  when 
significant  aliasing  began.  This  indicates  that  very  little  energy  leaked  off  the 
spatial  grid  due  to  beam  divergence.  Actually,  the  use  of  Fourier  transforms  in 
the  propagation  preserved  the  energy  so  that  any  energy  that  leaked  off  the 
grid  was  aliased  back  into  the  field  on  the  grid.  Energy  that  leaks  off  in  the 
spatial  domain  results  from  undersampiing  in  the  spatial  frequency  domain,  and 
vice  versa. 

Aliasing  most  often  occurs  from  undersampiing  in  the  spatial  domain, 
which  corresponds  to  energy  leaking  off  the  grid  in  the  spatial  frequency 
domain.  Prudent  choice  of  source  and  the  final  irradiance  patterns  across  the 
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grid  can  minimize  the  amount  of  energy  ttiat  leaks  off  the  grid  to  be  aliased 
back  in  for  the  spatial  domain.  However,  as  turbulence  strength  increases,  the 
coherence  length  decreases  and  more  energy  occurs  at  higher  and  higher 
spatial  frequencies.  Eventually  spatial  undersampling  occurs  and  this  high 
spatial  frequency  energy  leaks  off  the  spatial  frequency  grid  and  is  aliased  back 
in.  Figure  (22)  gives  the  radial  power  spectral  density  corresponding  to  Fig. 

(21)  and  shows  the  spread  of  energy  to  higher  spatial  wavenumbers  as 
turbulence  strength  increases.  The  low  turbulence  spectra  possess  a  strong 
central  peak  that  becomes  more  rounded  and  flatter  as  turbulence  strength 
increases,  causing  more  energy  to  leak  off  the  grid.  Even*  jally,  enough  energy 
has  leaked  off  the  grid  and  aliased  bade  in  to  make  the  spectrum  approximately 
uniform  at  all  spatial  frequencies. 

Figure  (23)  actually  shows  this  energy  being  reflected  back  in  after  it 
spills  off.  To  reveal  this  phenomenon,  the  radial  power  spectrum  for  a  512x512 
grid  has  been  divided  by  the  corresponding  part  of  the  power  spectrum  on  a 
1024x1024  grid.  Since  the  larger  grid  has  been  chosen  with  a  finer  mesh  than 
the  512x512  grid,  the  larger  grid  will  not  experience  significant  aliasing  as  soon. 
As  turbulence  strength  increases,  the  ratio  shows  the  enhancement  of  energy 
at  the  edge  of  the  512x512  grid  (bins  200-256)  as  the  energy  leaks  off  /  reflects 
back  in.  At  very  high  levels  of  turbulence,  the  aliased  energy  has  spread 
inward  over  all  spatial  frequencies. 
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Power 


FIgm  22  Power  spectrum  radial  profile  as  turbulence  strength  increases:  p, 
-  0  (solid  line),  0.5  (pluses),  1.5  (asterisks),  5  (dots),  15  (diamonds),  50 
(triangles),  150  (squares),  500  (X's). 


0  so  100  ^50  200  250 
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Figure  23  Power  spectrum  radial  profile  for  512x512  grid,  divided  by 
1024x1024  profile,  showing  energy  aliasing  back  in.  0.15(plus), 
0.5(asterisk),  1 .5(diamond),5(triangle).  1 5(square),50(X), 1 50(dot). 
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Aliasing  appears  to  causa  the  intensity  pe^ng  behavior  of  the  central 
field.  These  investigations  parameterized  this  onset  of  peaking  with  five 
methods,  which  provided  guidelines  for  maximum  turbulence  strength 
valid  for  a  given  grid  size.  An  irradiance  ratio  and  a  Strehl  ratio  were  first  used 
to  identify  the  onset  of  the  peaking.  A  Gaussian  source  was  propagated  at 
turbulence  strengths  ptt*  =  [5x10^.  5x10-®.  5x10-*,  0.15,  0.5,  1.5,  5,  15,  50, 
150,  500]  using  64x64,  128x128,  256x256,  512x512,  and  1024x1024  grids. 

Ten  runs  for  each  case  were  used  to  caiculate  the  average  irradiance  as  a 
function  of  radius.  These  irradiance  profiies  were  used  to  calculate  an 
irradiance  ratio,  defined  as  the  irradiance  at  the  grid  center  divided  by  the 
irradiance  at  the  maximum  grid  radius,  and  also  the  Strehl  ratio,  defined  as  the 
irradiance  at  the  grid  center  for  the  E-field  propagated  through  turbulence 
divided  by  the  center  irradiance  for  an  E-fieid  propagated  through  zero 
turbuience.  Figure  (24)  plots  the  irradiance  ratio  versus  turbulence  strength  for 
different  grid  sizes,  and  Fig.  (25)  shows  the  Strehl  ratio  versus  turbulence 
strength  for  different  grid  sizes. 

As  turbulence  strength  increased,  diffraction  and  scattering  spread  the 
energy  outward  and  caused  the  irradiance  and  Strehl  ratios  to  decrease. 
Eventually,  both  ratios  reached  a  minimum  and  then  started  increasing  because 
of  the  alias-induced  intensity  peaking.  This  minimum  occurred  at  higher 
turbulence  strengths  for  larger  grid  size  because  the  larger  grid  sizes  sampled 
with  a  finer  mesh  and  did  not  significantly  alias  as  soon.  Using  estimates  of 
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Figura  24  Irradiance  ratio  versus  turbulence  strength  (  =  [5x1  5x10^ 

for  grid  sizes;  64x64  (asterisks),  128x128  (diamonds),  256x256  (triangles), 
512x512  (squares),  1024x1024  (X's). 
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these  minima  as  a  guide  to  the  onset  of  significant  aliasing,  Fig.  (26)  plots 
these  Po’  values  versus  grid  size  from  the  irradiance  and  Strehl  ratios,  along 
with  least  squares  fits  and  extrapoiations  to  larger  grid  sizes. 

Martin  and  Fiattd  (1988)  and  others  have  validated  computer  simulations 
for  predicting  statistical  properties  such  as  normalized  irradiance  variance  from 
the  E-fieid  propagated  through  turbulence.  The  departure  of  the  computer 
simulation  E-fields  and  intensities  from  their  known  or  expected  smooth 
behavior  as  turbulence  increases  represents  another  telltale  sign  of  significant 
aliasing.  Figure  (27)  shows  the  normalized  irradiance  variance  calculated  from 
computer  simulated  E-fields  versus  po’  using  grids  of  sizes  64x64,  128x128, 
256x256,  512x512,  and  1024x1024.  The  irradiance  value  used  for 
normalization  was  taken  as  the  average  irradiance  over  the  central  calculation 
region.  Though  this  single  value  normalization  is  not  optimal  (see  the  following 
section,  F.  AddWonal  Simulation  Parameters),  it  does  make  the  normalized 
irradiance  variance  calculation  sensitive  to  the  peaking  behavior  because 
peaking  adds  large  amounts  of  unphysical  variance.  In  the  Rytov  regime 
(3o^  ^  1.0),  all  grids  successfully  simulated  the  E-field  as  indicated  by  their 
normalized  irradiance  variances  agreeing  with  the  Rytov-Tatarski  theory.  In  the 
saturation  regime  (Pg^  ^  1.0)  the  normalized  irradiance  variance  saturates  and 
turns  downward,  eventually  approaching  unity  according  to  asymptotic  theory. 
Assuming  that  the  1024x1024  grid  provides  the  most  accurate  propagation 
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Figure  27  Normalized  inradiance  variance  versus  turbulence  strength  from 
simulations  with  different  grid  sizes. 
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simulation,  the  64x64  grid  normalized  irradiance  variance  departs  from  the 
1024x1024  values  before  the  peak  was  reached,  and  128x128  grid  normalized 
irradiance  variance  departed  just  beyond  the  peak.  Their  meshes  were  not 
small  enough  to  adequately  sample  the  E-field,  and  the  resulting  aliasing 
caused  peaking  that  drove  up  the  normalized  irradiance  variance  calculation. 
The  256x256  grid  produced  the  saturation  peak,  but  soon  suffered  from 
significant  aliasing.  The  512x512  calculations  successfully  produced  the  peak 
but  departed  from  the  1024x1024  predictions  before  Po*  ”  50  .  Judging  from 
these  behaviors,  the  1024x1024  grid  probably  remains  valid  through  3,,^  >  50  . 
All  grids  eventually  showed  the  normalized  irradiance  variance  anomalously 
rising  for  high  enough  turbulence  strength  because  aliasing  produced  peaking. 
Estimates  of  the  maximum  values  at  which  the  five  grid  sizes  remained  valid 
are  plotted,  along  with  a  least  squares  fit  extrapolated  to  larger  grid  sizes,  in 
Fig.  (28). 

In  a  similar  manner,  the  coherence  length  and  the  half  width  at  half 
maximum  (HWHM)  of  the  atmospheric  MTF  predicted  using  simulations  with 
different  grid  sizes  also  show  unphysical  behavior  as  aliasing  became 
significant.  Figure  (29)  plots  the  simulation  derived  coherence  length, 
normalized  by  the  Kolmogorov  turbulence  theoretical  values,  versus  turbulence 
strength  measured  by  p^^.  Figure  (30)  shows  the  corresponding  plot  for 
HWHM.  Both  show  the  eventual  strong  rise  in  coherence  length  and  HWHM 
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FigiNV  29  Spherical  wave  coherence  length  from  simulations  versus  turbulence 
strength  for  grid  sizes  128x128  (diamonds),  256x256  (triangles),  512x512 
(squares),  and  1024x1024  (X's). 
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Figure  30  Spherical  wave  HWHM  from  simulations  versus  turbulence  strength 
for  grid  sizes  128x128  (diamonds).  256x256  (triangles).  512x512  (squares),  and 
1024x1024  (X’s). 
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compared  to  theory  for  strong  turbulence  when  aliasing  produces  peaking. 
Estimates  of  the  minima  of  the  plots  can  be  used  as  another  indicator  of  the 
turbulence  strength  at  which  significant  aliasing  occurs.  Figure  (31)  plots  these 
points  versus  grid  size,  along  with  a  least  squares  estimate  extended  to  larger 
grid  sizes. 

The  onset  of  peaking  and  departure  from  expected  physical  behavior 
provide  symptoms  of  aliasing  but  do  not  indicate  the  amount  of  aliasing 
required  to  cause  them.  Determination  of  the  fraction  of  total  field  energy  that 
is  aliased  serves  as  one  parameterization  of  the  amount  of  aliasing  and 
accomplishes  two  goals:  (1)  relate  the  amount  of  aliasing  occurring  to  an  easiiy 
measurable  characteristic  of  the  simulation  ,  and  (2)  determine  how  much 
aliasing  must  occur  to  invalidate  ttie  computer  simulation. 

First,  to  parameterize  the  amount  of  energy  aliased,  the  same  size 
Gaussian  source  was  applied  to  64;^,  128x128,  256x256,  and  1024x1024 
grids  and  then  propagated  through  turbulence.  The  width  of  the  source  was 
chosen  so  that  the  Fourier  transform  had  approximately  the  same  width  on  the 
final  grid  and  the  Gaussian  sources  on  the  different  grids  were  normalized  to 
the  same  energies.  Since  identical  sources  were  applied  in  the  spatial  domain 
according  to  the  analytical  Gaussian  formula,  the  spatial  frequency 
representation  of  these  sources  improved  as  the  grid  size  increased  and  the 
mesh  became  finer.  The  1024x1024  grid  served  as  an  approximation  to  infinite 
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Figure  31  Estimates  and  least  squares  fits  from  spherical  wave  coherence 
length  (triangles  /  asterisks)  and  HWHM  (diamonds  /  pluses)  of  the  maximum 
turbulence  strength  versus  grid  size  for  valid  simulations. 


grid  size.  Due  to  its  larger  extent  and  finer  mesh,  the  1024x1024  spatial 
frequency  representation  of  the  Gaussians  contained  spectral  energy  in  the 
region  outside  the  spectral  footprints  of  the  smaller  grids,  as  Fig.  (32)  illustrates. 
Since  each  grid  started  with  the  same  total  energy,  the  fraction  of  spectral 
energy  in  the  1024x1024  grid  lying  outside  the  spectra)  footprints  of  the  smaller 
grids  approximated  the  amount  of  energy  aliased  in  the  smaller  grids.  This 
process  was  then  applied  to  propagations  of  the  Gaussian  beams  through 
increasing  strengths  of  turbulence,  well  into  the  significant  aliasing  regions  for 


Figtae  32  Illustration  of  the  spectral  representations  of  the  Gaussian  beam 
with  different  grid  sizes. 
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the  smaller  grids.  This  method  provided  estimates  of  the  fraction  of  energy 
aliased  versus  turbulence  strength  for  the  three  smaller  grids. 

To  achieve  the  first  goal  of  relating  the  amount  of  aliasing  to  a 
measurable  simulation  parameter,  the  average  radial  power  spectral  density  of 
the  propagated  E-field  was  calculated  for  each  turbulence  strength  with  the 
smaller  grids.  Using  these  profiles,  a  spectral  ratio  was  calculated,  defined  as 
the  ratio  of  the  power  spectral  density  at  the  grid  center  to  the  power  spectral 
density  at  the  maximum  grid  radius.  Figure  (33)  plots  these  results,  and  also 
includes  spectral  ratios  for  512x512  and  1024x1024  grids.  Using  the 
information  on  fraction  of  energy  aliased  versus  and  spectral  ratio  versus 
Po^.  the  turbulence  strength  p,,^  served  as  the  connection  between  fraction  of 
energy  aliased  and  the  measurable  simulation  parameter  of  spectral  ratio. 

Figure  (34)  shows  plots  of  these  spectral  ratios  versus  corresponding  fraction  of 
energy  aliased.  The  log-log  plot  behavior  proved  approximately  linear  over  the 
range  0.001  to  0.1  ,  which  appears  very  useful  because  0.1%  of  energy  aliased 
probably  does  not  affect  the  simulation  results  significantly  while  more  than 
10%  of  energy  aliased  probably  does. 

Using  the  data  points  of  Fig.  (34)  and  performing  a  linear  least  squares 
fit  to  the  logarithms 

loatoCfl) '  eiog«l«5)  *  iobwA  (ii«) 
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or  equivalently, 


/?  -  10^  O'***) 

where  R  refers  to  spectral  ratio  and  FR  refers  to  fraction  of  energy  aliased, 
gives  A-  -1.1  ±  0.1,  B-  -2.1  ±  0.1  .  Substituting  the  A  and  B  values  into  Eq. 
(Ill)  and  rearranging. 


1 


(112) 


This  equation  achieves  the  first  goal  of  relating  the  fraction  of  energy  aliased  to 
a  measurable  quantity  from  the  simulation. 

However,  Eq.  (112)  should  be  used  with  caution  because  it  appears  to 
be  specific  to  the  narrow  Gaussian  source  used  to  derive  it.  Figure  (35)  shows 
the  same  calculations  carried  out  for  a  source  that  was  the  Fourier  transform  of 
an  aperture  of  radius  equal  to  3/4  of  the  64x64  grid  radius  (i.e.  Airy-type 
source).  The  linear  region  did  not  appear,  though  the  plot  matched  the 
Gaussian  source  plot  when  the  fraction  of  energy  aliased  was  greater  than 
0.1  .  The  differences  arise  from  the  different  spectral  representations  of  the 
propagated  E-fields.  The  abrupt,  cylinder-shaped  irradiance  patterns  start  with 
more  energy  at  higher  spatial  frequencies  in  the  spectral  domain  than  the 
Gaussian  patterns  and  thus  never  get  below  0.001  fraction  of  energy  aliased 
even  at  very  low  turbulence  strengths  and  high  power  spectral  ratios. 
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Flgim  35  Power  spectral  density  ratio  versus  fraction  of  energy  aliased  for 
Airy-type  source;  grid  sizes;  64x64  (dotted  line),  128x128  (dashed  line), 
256x256  (solid  line). 
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The  sea>nd  goal  of  determining  the  amount  of  aliasing  that  invalidates 
the  simulation  used  the  fraction  of  spectral  energy  aliased  to  generate  another 
guideline  for  maximum  strength  of  turbulence  for  a  given  grid  size  where  the 
maximum  fraction  of  energy  allowed  to  be  aliased  in  a  simulation  was  10%. 
Equation  (112)  predicts  a  spectral  ratio  of  approximately  10  for  10%  of  energy 
aliased.  Linearly  interpolation  between  the  data  points  of  Fig.  (33)  provided 
estimates  of  the  turbulence  strengfris  C„^  ( thus  3^^ )  that  gave  a  spectral  ratio 
=  10,  for  10%  energy  aliased,  for  the  five  grid  sizes.  Figure  (36)  plots  these 
values  and  a  least  squares  fit  extended  to  larger  grid  sizes.  The  least  squares 
fit  as  a  function  of  grid  size  N  was 

•ogioPo  =  0.9  log,o(  V )  -  0.9.  ('•‘*3) 

or, 

p5  « 0.1 

The  fraction  of  energy  aliased  also  provided  estimates  of  maximum 
for  the  Airy-type  source.  Figure  (37)  plots  the  fraction  of  energy  aliased  versus 
turbulence  strength  for  the  64x64,  128x128,  and  256x256  grids.  Again 
choosing  10%  energy  aliased  and  using  linear  interpolation  yields  data  points 
for  the  solid  line  in  Fig.  (36)  that  shows  the  maximum  turbulence  strengths 
obtained  from  Fig.  (37).  For  10%  energy  aliased,  the  predictions  for  maximum 
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Figure  36  Least  squares  fit  (pluses,  dotted  line)  of  maximum  beta  squared 
versus  grid  size  for  a  Gaussian  source  that  had  10%  energy  aliased  and  for  an 
Airy-type  source  (asterisks,  solid  line). 


89 


0.000 1 


0.0010  aoloo  0.1000  i.oooo  lo.oooo 

0«to  squorad 


100.0000 


npM  37  Frtcoon  of  onorgy  linod  vorous  turtufonoo  ttrongth  wtth  an  Airy 
typa  touroa  for  grid  tiZM  64x64.  121^126.  and  256x256 


for  the  three  smeller  grids  using  the  Airy-type  source  agree  within  ~'20%  with 
the  predictions  fror«i  the  Gaussian  source.  Lower  values  of  frac^onal  energy 
aliased  do  not  provide  such  agreement. 

Another  measure  of  maximum  turbulence  strength  involved  the 
coherence  length  and  grid  size.  Aliasing  occurs  because  the  E-fieid  fluctuates 
significantly  over  scale  sizes  smaller  than  the  grid  element  size.  Intuitively, 
some  conriection  should  exist  between  coherence  lerHJth  of  the  E-field  and  foe 
onset  of  significant  aliasfog.  Fried's  coherence  length  r,  represents  the  distance 
over  which  the  atmospheric  MTF  falls  to  exp(-3.44)  »  0.032  (Fried,  1966,  p. 
1360-1383).  Equation  (38),  which  relates  to  C.^  and  Eq.  (66),  which 
relates  r,  to  C.^  arxl  Eq  (92),  which  relates  the  element  size  Ax  to  N, 
when  combined,  /ield  foe  turbulence  tinngHh  correspondir>g  to  a  given  grid  size 
for  a  specific  number  (y)  of  Ax's  per  r, 

.O.m  ,9pfm1ot/m¥9 

pj  •  I 

'  0429  m¥9 

Figure  (38)  plots  these  p,’  for  y  «  2.  2.5,  and  3  for  foe  sphericai  wave  case 
The  10%  aliaaed  energy  Mne,  indicated  by  plus  symbols,  lies  newest  foe  2.5  Ax 
per  r,  Mne  The  E-fieW  must  be  spabaiy  swT^)led  with  2.5  Ax  per  r,  to  limit  foe 
flracbon  of  energy  aNased  to  10%,  analogous  to  foe  Nyquist  criforion  two 
samplM  per  cycle 
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(39)  provides  a  combined  piot  of  tt)e  maximum  turbulence 
strength  that  produced  valid  E-fields  for  a  given  grid  size  from  all  of  the 
fore^ing  measures  isity  ratio.  Strehl  ratio,  normalized  irradiance  variance, 
coherence  length,  HWhM,  fraction  of  energy  aliased,  and  coherence  length 
(y  s  2.5).  Most  estimates  tie  witNn  a  factor  of  2  of  each  other,  and  the  10% 
energy  aliased  line  (dotted/pluses)  given  by  Eq.  (114)  represents  an 
approximate  lower  bound  to  those  estimates  based  upon  onset  of  significant 
aliasing. 

Figure  (39)  provided  three  s«,|r^.  .4  corKlusions:  (1)  the  computer 
simulations  remained  v^  until  approximatei}  10%  of  the  energy  became 
aliased  (achieving  the  ffrst  goal),  (2)  using  the  10%  energy  attased  Kne,  the 
maximum  turbulence  strength  p,’  for  valid  E-flelds  for  s  grid  size  M  was 

and  (3)  maxtonum  valid  turbuierKe  sfrength  corresporxfed  to  approximately 
2.5  grid  elements  per  r,  (based  on  prox^ity  of  the  y  -  2.5  line  to  the  10% 
energy  aliased  Hne).  CocKlusion  (2)  implies  that  doubto>g  the  grid  size  N  (with 
grid  element  size  Ax  given  by  Eq.  (92))  sUghtty  less  than  doubles  the  maximum 
Pq’  that  the  can  simirfate,  and  that  these  investigations,  which  use  a 
1024x1024  grid,  shO(4d  be  vaM  up  to  p,*^  50. 


93 


F.  ADDITIONAL  SIMULATION  PARAMETERS 

Additional  aspects  of  a  computer  simulation  were  addressed  to  ensure 
validity,  including  the  number  of  phase  screens,  the  method  of  normalizing  the 
irrat^ance  variance,  the  number  of  realizatkms  (propagations  through 
turbulence)  required  for  representative  statisbcs,  arKi  the  width  of  the  final 
irra<fiance  field  on  the  grid. 

The  number  of  phase  screens  must  be  large  enough  to  represent  the 
turbulence  accurately  along  the  path  and  produce  proper  irradiance  and 
coherence  statistics.  Martin  and  FlatM  (1988)  determ^ied  the  number  of  phase 
screens  by  requiririg  that  the  variance  due  to  propagation  over  the  distance  Az 
between  phase  screens  be  less  than  1/10  the  vwianoe  from  the  total 
propagation  over  the  distance  L 

of(Az)  <  ai  oJ(L),  (1^7) 

and  additionally  that  the  vMue  of  the  variance  from  one  step  be  less  than  0. 1 

o^{Ai)  <  0.1  . 

WHh  these  coruideratlorui.  they  generally  used  20  phase  screeiis  for  the^ 
simulations 

These  investigations  examir>ed  the  irradianoe  and  coherence  statistics 
(firecfiy  arKi  concluded  that  approximately  30  phase  screens  were  recMr*0  to 
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ensure  simulation  validity.  Figure  (40)  shows  a  set  of  spherical  wave 
normalized  irradiance  variance  simulations  for  a  512x512  grid,  ~  1.5,  and 
the  number  of  phase  screens  varying  between  2,  4,  8,  16,  32,  and  64. 

Assuming  the  64  phase  screen  case  most  doseiy  approximated  physical  reality, 
as  UfH  as  8  phase  screens  gave  a  normalized  irradiance  variance  within  ~10%  of 
this  value.  By  32  phase  screens,  the  ncxmalized  irradiance  variance  had 
stabilized  to  within  2%  of  the  64  phase  screen  value.  For  this  level  of 
titfbulence  that  lies  at  the  beginnirtg  of  the  saturation  regime,  the  addition  of 
phase  screens  beyond  32  did  not  affect  the  normalized  irradiance  variance 
significantly. 

Figure  (41)  shows  the  normalized  ^radiance  variance  from  propagations 
with  2,  4,  8,  16,  32,  and  64  phase  screens  at  50  with  a  1024x1024  grid 
and  30  realizations  for  each  manber  of  phase  screens.  This  strength  of 
turbulence  represents  the  limit  of  simulation  validity  with  the  1024x1024  grid 
and  lies  in  the  strong  saturation  regime  beyond  the  normalized  irradiance 
variance  peak.  In  this  case,  32  phase  screens  provided  a  normalized 
irradiance  variance  within  3%  of  the  64  phase  screen  case,  arKl  the  16  phase 
screen  value  was  within  approximately  10%.  Because  32  phase  screens 
appears  to  increase  the  accuracy  of  the  simutation  by  ~5%  compared  to  16 
phase  screens,  me  larger  number  of  phase  screens  was  chosen  for  these 
investigabons 
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Plgum  41  NomwIizoO  irradlanot  vartano*  as  a  function  of  numOar  of  phase 
screens  for  turbulence  strength  p,’ «  SO. 
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Using  th«  sssne  1024x1024  reslizattons,  the  coherence  length  and 
HWHM  of  the  atmospheric  MTF  were  calculated  for  the  2,  4,  8,  16,  32,  and  64 
phase  screen  cases.  Figure  (42)  plots  the  results.  For  coherence  length,  the 
32  phase  screen  value  was  within  approximately  1%  of  the  64  phase  screen 
value,  and  the  16  phase  screen  value  was  within  approximately  3%.  The 
HWHM  plot  indicates  2%  ar>d  3%  agreements  for  32  and  16  phase  screens, 
respectively.  The  32  phase  screens  used  for  these  investigations  thus  proved 
sufncient  for  cohererKe  length  and  HWHM  calculations. 

The  method  of  normaltzing  the  irra<tonce  variance  and  the  number  of 
realizations  to  use  for  statistical  accuracy  proved  to  be  interrelated 
considerations  Turbulence  (ftffracts  and  scatters  energy  outward  from  an 
initialty  well-defined  beam  For  the  Axy-type  source  used  in  these 
investigations,  the  average  irradiance  over  the  central  portion  of  the  final 
propagated  field  was  uniform  for  zero  turbulence  but  became  a  function  of 
radial  cfistance  from  the  propagation  axis  when  turbi^erKe  was  present.  This 
was  an  artifact  of  foe  computer  simulation  that  had  to  use  a  beam  spatially 
cortfined  to  the  grid  rather  than  a  true  spherical  (or  plane)  wave  Figure  (43) 
shows  irradiance  versus  radial  distance  r  for  foe  ceritral  256x256  portion  of 
1024x1(^4  (fod  simulations,  averaged  over  30  reaitzations,  and  using  a 
turbulence  strength  «  50  The  average  radM  irradlWKe  vvies  by 
iHPproximafoly  15%  ovw’  the  widfo  of  this  calculation  region  This  rarfial 
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variation  of  the  average  irradiance  was  removed  from  the  rKKmaltzed  irradiance 
variance  calculations  with  a  method  smilar  to  that  of  Flattd,  Wang,  and  Martin 
(1993).  The  calculation  region  was  restricted  to  half  the  radius  of  the  zero 
turbulence  irradiance  pattern,  which  was  previously  chosen  as  one  half  the  grid 
radius.  WHh  the  1024x1024  grid,  this  calculation  region  equaled  the  largest 
radius  circle  inside  the  central  256x256  portion  of  the  grid.  This  disk  was 
further  divided  into  concentric  rings  arKl  the  irradiance  in  each  ring  averaged 
over  ail  30  realizations.  These  average  ring  irradiar)oes  were  used  to  normalize 
the  irradiance  variar>ce  calculated  from  each  field  Smaller  ring  size  (thus  more 
rings)  reduced  the  number  of  poirrts  in  the  irradianoe  average  for  each  nng  and 
required  more  realizations  to  ensure  a  sufficient  nunber  of  points  to  yield  a 
stable  average  irradiance. 

The  rir>g  width  had  to  be  smaM  enough  to  compensate  for  the  radial 
variation  in  average  irradianoe  and  the  number  of  realtzations  targe  enough  to 
provide  enough  points  to  yield  represerrtative  average  values  To  determine 
suitable  ring  width  and  number  of  realizations.  50  1024x1024  realizations  with 
the  Airy-type  source  at  »  1.5,  3,  10.  arKt  20  were  run  and  the  normakzed 
irradiance  variances  caloMted  usirtg  ring  widths  of  128,  32.  8,  4,  2,  arxl  1  grid 
elements,  dx,  and  with  number  of  ra^zations  in  the  average  equal  to  1,  5,  10, 
15,  .  ,50  Figure  (44)  plots  the  rest^  for  the  89*  s  10  set  The  pluses  lir>e 
corresponds  to  1  ring  of  width  128  dx  (i  e  the  whole  calculation  region)  and  a 
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8tabi«  rvormaltzed  irr»tianc8  value  was  achieved  using  about  20  realizations. 
However,  a  single  hrtg  does  not  compensate  for  the  radial  variation  in  average 
irra<^ar)ce  and  can  thus  give  erroneous  normalized  irradianoe  values  (for  the 
same  reason,  it  becomes  sensitive  to  the  pe^ng  behavior  from  significant 
aliasing)  The  X  line  represents  foe  division  of  foe  disk  into  4  rings  of  width  32 
Ax  and  recM^  ‘‘25  realizations  in  foe  average  to  give  a  stable  normalized 
irradiance  vwiance.  SmaHer  ring  widths  (Unas  with  squares,  diamonds, 
triangles,  and  asterisks  lines)  all  showed  similar  behavior  and  required  -  30  runs 
to  reach  a  stable  average  normi^ed  irradiance  variance.  The  =  1.5,  3, 
arrd  20  runs  aN  provided  similar  results.  A  similar  p,*  s  io  series  with  a 
Gaussian  source  required  a  mir^mum  of  8  rings  (16  Ax  each)  and  30 
re^izations  to  active  stable  average  normalized  irradianoe  variarK^es. 
Consec^jentty,  these  investigations  used  30  raakzations  as  a  guideline  for  valid 
simi4ation,  and  chose  a  rirrg  width  elements  to  irilow  foe  greatest 

adjustment  to  real  variations  in  average  radial  foterrsity  and  yet  remain 
computationally  efficient.  As  Fig.  (44)  ^lows.  insufficient  averaj^  over  enough 
realizations  to  yield  truly  representative  average  values  can  reduce  foe 
normalized  irrackanoe  varUmca  by  15  •  30%. 

Sknkarty,  stable  coherence  length  values  (discussed  in  section.  H. 
Cebetetice  Langfo)  required  averaging  over  multiple  realizations.  To  determine 
the  number  of  realizations  requked  fo  foe  average,  50  realizations  usfog  a 
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1024x1024  grid  were  g«n«r«t«d  for  p,’  s  3  and  20  and  the  central  256x256 
portion  of  the  fields  again  used  to  calculate  the  atmospheric  MTFs.  These 

individual  MTFs  were  then  averaged  in  groups  of  1,  5,  10,  15 . 50  before 

calculating  the  coherence  ler>gths.  Figure  (45)  shows  the  resulting  coherence 
lengths,  r,.  (These  particular  runs  used  L  >  150  m  and  X  =  420  nm,  giving  a 
FresTHil  length  ■  L  /  1  %  -  10  mm  and  resultir>g  in  millimeter-sized 
colwence  lengths.)  As  few  as  5  realizations  in  the  average  yielded  coherence 
lengths  within  4%  of  the  50-realization  values.  These  investigations  still  used 
30  realizatiorM  in  the  average  because  these  30  fields  were  already  available 
from  the  normalized  irradiance  variance  calculations. 

The  irradiarKe  ratio  provided  an  indicator  of  the  appropriate  final  beam 
radius  to  use  to  avoid  excessive  scatter  of  enwgy  off  the  grid.  The  maximum 
Pq’  for  a  given  grid  varies  with  the  final  beam  radius  because  larger  radii  scatter 
energy  off  the  grid  sooner  and  cause  aliasmg  at  lower  To  characterize  this 
behavior,  64x64,  128x128,  and  256x256  grid  simulations  were  run  at  turbulence 
*<>’*hgths  (5x10'*,  5x10^  for  ffoal  bewn  radii  of  4/8,  5/8,  6/8,  and  7/8  grid 
radius,  and  the  resulting  irradiartce  ratios  (average  irradiance  at  center  divided 
by  average  irracttarKe  at  grid  radius)  were  plotted  versus  ^  .  Figure  (46) 
shows  the  plot  for  the  256x256  grid.  Again,  the  mir)ima  correspond  to  the  onset 
of  peaking  and  occur  at  lower  Pp’  for  larger  final  beam  radius.  Linear  least 
squares  fit  of  the  final  beam  radii  to  these  V  values  for  the  minimum 
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irradianoe  ratios  indicates  that  a  final  beam  radius  of  approximately  0.7  grid 
radius  corresponds  with  the  10%  energy  aliased  cutoff  for  the  grids.  A  final 
beam  of  0.7  grid  radius  provided  the  largest  illuminated  central  region  while  sfill 
meeting  the  10%  aliased  energy  criterion.  To  be  somewhat  conservative,  these 
investigations  used  a  final  beam  radius  =  0.5  grid  radius  to  reduce  further  the 
energy  scattered  off  the  grid. 

Q.  PHASE  SCREENS 

1.  Pliaae  Saaen  Oeneratioii 

WHh  the  split-step  method,  the  effects  of  turbulence  along  the 
optical  path  are  introduced  into  the  simulation  by  dividing  the  optical  path  into 
steps  and  applying  a  rarKfom  phase  to  the  complex  E-fieid  at  each  step.  As 
long  as  the  steps  are  smalt  erKXjgh  that  geometrical  optics  approximately 
applies,  the  E-fleld  only  acquires  a  random  phase  change  as  it  propagates 
across  each  step  (Knepp,  1983).  Diffraction  as  the  field  propagates  across 
many  steps  then  produces  the  ampfitude  variations.  The  random  phases  are 
assumed  to  be  Gaussian  distributed  about  a  zero  mean  with  variance 
proportionai  to  the  turbulenoe  strength  C,’  and  possessing  spatial  structure 
function  consistent  with  the  assumed  Kolmogorov  turbulence  and  appropriate 
inner  scale.  Knepp  (1963)  and  Martin  and  Flatti  (1988)  describe  the  process 
of  generating  the  phase  screen  with  these  characteristics,  as  discussed  below. 
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The  phase  screen  generation  begins  in  the  spatial  frequency 
domain  by  imposing  the  proper  spatial  structure  function.  An  NxN  grid  of 


complex  numbers  6o(K„Ky)  is  formed  whose  real  and  imaginary  parts  are  each 
Gaussian  distributed  random  numbers  with  zero  mean  and  unity  standard 
deviation.  This  6o(K,.Ky)  represents  the  Fourier  transform  of  a  grid  of 
uncorrelated  Gaussian  distributed  random  numbers  6o(x.y)  representing  phases. 
The  proper  spatial  structure  function  corresponding  to  turbulence  statistics  is 
imposed  upon  the  random  phases  6o(x.y)  by  applying  a  filter  A(K,.iCy)  to 

e(lC,.K,)  »  A{Kg,Ky)  (Il#) 

Taking  the  magnitude  of  both  sides  of  Eq.  (119),  squaring,  assuming  that  the 
filter  function  is  real,  and  then  takirtg  expectation  values  gives 

<  I  e(v,Ky)  I*  >  •  <  I  I*  >  • 

Where  use  was  made  of  the  fact  that  the  (3aussian  random  numbers  6o(K„Ky) 

have  a  variance  of  1.  The  two-dimensional  power  spectral  density  of  the  phase 

^s(k»S)  rela^  to  0(k„k,)  by  (Goodman,  1985) 

Where  Ak  represents  the  grid  element  size  in  the  spatial  frequency  domain. 

The  Hankei  transform  of  the  power  spectra  density  F,(k)  gives  the  phase 
structure  function  D,(p)  that  characterizes  the  spatial  distribution  of  the  phase 
fluctuations  of  the  E-field  (Tatarski,  1961) 
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0»(P)  •  /1 1  -  4(*P)  1  ''»(*.0)  K  <fc. 

-«•  _ 

Where  local  iaotropy  has  been  assumed.  Tatarski  derived  the  relation  between 
the  two-dimensional  power  spectral  density  of  phase  fluctuations  F,(ic„Ky)  wkI 
the  tiwee-dimensional  spectrum  of  index  of  refraction  fluctuation-.  <b,(K), 

-  2»lf*  L  ♦^(k  ),  (123) 

where  x  •  ♦  kJ  .  This  sequence  of  steps  means  that  the  phase  screen 

can  have  the  proper  spatial  statistics  by  starting  with  the  proper  spectrum  of 
refractive  index  fluctuations  (hence,  the  proper  structure  function). 

The  spectrum  of  index  of  refraction  fluctuations  assuming 
Kolmogorov  turbuienoe  with  inner  scM  is 

-  OJOBi  C,*W  «■”"  P(ieW.  (1W) 

Where  F(k^  gives  the  inner  scale  d<^>endence  for  kmer  scale  ^  (see  Fig.  8). 
Siriistituting  Eq.  (124)  into  Eq.  (123)  gives  the  power  spectral  density  of  phase 
fluctuations 

Fg{x^)  -  2iclr*  L  (0.033)  Ci{z)  F(k^),  (««) 

and  using  Eq.  (121)  spedfles  the  proper  form  for  6(k,.k,),  the  corresponding 
Fourier  transform  of  the  random  phases 

<  I*  >  *  (As)-*  2xk*  L  (04»3)  C!{J^  F(k^).  (12«) 


Equation  (120)  than  gives  the  corresponding  filter  function  A(K.,Ky)  to  apply  to 
the  array  of  complex  numbers  6o(ic„tCy) 

A(ic„Ky)  »  (Ax)-^  L  (0.033)  C!!{z)  F(ic^). 

Since  the  Kolmof^rov  spectrum  (x  has  a  singularity  at  k  -  0,  ttie 
point  in  the  filter  fur>ction  is  set  to  zero,  removing  overall  piston  (i.e.  common 
phase  offset  over  whole  screen)  from  the  phase  (Cochrane,  1985)  and  keeping 
the  spectral  energy  finite.  Conversion  to  a  discrete  grid  representation  occurs 
by  substituting  x,  =  n.  Ax,  k,  =  n^  Ak.  x*  *  (Ax)*  (n,*  ♦  n^*),  arnl 
Ax  2ic/(N  Ax),  where  Ax  is  the  grid  element  size  in  the  x-y  domain  and  (n,.ny) 
are  grid  coordinates.  The  Fourier  transform  (FT)  of  the  filtered  array  of  random 
variables  gives  the  phase  screen  in  the  spatial  domain 

e(x,y)  .  0.0864  Hr  L  {N  Ax)^  FT[  ]• 

Since  this  phase  screen  is  actually  complex-valued,  both  the  real  and  imaginary 
parts  represent  vaNd  random  phase  screens  that  were  separately  applied  to  the 
E-field  (Cochrane.  1965). 

2.  Low  8pa6al  Frequency  Conecdon 

Due  to  the  finite  size  of  the  the  above  phase  screen  will  not 
have  the  proper  structure  function  for  separations  of  the  order  of  the  grid  width. 
Low  spatiai  fTsquency  components,  espedatiy  tilt-type  terms,  are  under- 
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represented  (Cochrane,  1985).  These  computer  simulations  incorporate  an 
algorithm  formulated  by  Cochrane  that  employs  an  expansion  of  the  phase 
screen  in  a  Karhunen-Loeve  basis  set  (whose  components  are  mutually 
orthogonal)  to  correct  some  of  the  low  spatial  frequency  terms. 

The  general  idea  of  the  low  frequency  correction  builds  upon  the 
above  procedure  to  generate  a  phase  screen.  The  low  spatial  frequency 
contribution  to  the  phase  screen  is  a  superposition  of  orthogonal  low  spatial 
frequency  terms,  just  like  the  superposition  of  complex  exponential  terms  by 
Fourier  transform  in  the  phase  screen  generation  process  above.  The  strength 
of  each  low  spatial  frequer>cy  term  is  a  Gaussian  random  variable  with  an 
appropriate  variance,  just  as  the  Gaussian  raryJom  numbers  above  were  filtered 
to  give  the  proper  varUmce  and  thus  determine  the  strength  of  the 
corresporKJing  complex  exponential  in  the  spatial  domain.  Thus,  the  two 
objectives  involve  finding  appropriate  set  of  orthogor^al  low  spatial  frequency 
functions,  and  determkiing  the  corresporKiing  variances  applinble  to 
atmospheric  turbulence. 

An  arbitrary  ftmction  can  be  expanded  in  terms  of  Karhunen- 
Loeve  hjnctions  that  are  orthogonal  by  definition.  To  determkie  a  set  of 
Karhunen-Loeve  functions  appropriate  for  Kolmogorov  turbulence,  Cochrane 
builds  on  the  work  of  Non  arKi  (1975)  considers  expansion  of  an  arbitrary 
function  ^r,0)  over  a  drcuiar  aperture  of  radius  R  in  terms  of  Zemike 
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polynomialt  (Bom  and  WoHd,  1970) 


♦(f.e)  -  Sa.z(p.e). 

whore  p  represents  the  normalized  distance  r/R,  a,  are  the  expansion 
ooefftcients,  and  Z,  are  the  Zemike  polynomials.  With  W(r/R)  representing  the 
aperture  function,  the  coefficients  tn 

•l  ■  fd*'  m'lft)  ♦(A») 

Noll  assumed  that  these  coefficients  were  (Gaussian  random  variables  with  zero 
mean  arxl  with  a  covariance 

( afaf)  -  /«!»  fdf'  tn())Wtp')  ^(p.»)  <  ♦(«?)  M'fp')  > 

Fourier  transforming  to  the  spatial  frequency  domain 

<<•/  •//«<«'  o;w  •(WA.ic'/R)  ciKO. 

Where  Q^k)  represents  the  Fourier  transform  of  the  |th  Zemike  polynomial,  and 
^k/R,k'/R)  represents  the  Kolmogorov  spectrum  of  phase  fluctuations.  Non 
analyticaily  performed  the  integrals  to  give  a  covariance  matrix  in  which  the 
terms  represent  the  expected  covariances  due  to  Kolmogorov  turtulerioe. 

Cochrane  (1985)  notes  that  the  Zemike  polynomials  cannot  be 
used  to  form  an  orthogonal  expansion  of  the  turbuteiKeKlMorted  phase 
because  the  expansion  coefficients  are  correlated,  indicated  by  nonzero  off- 
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diagonal  elements  In  Noll's  covariance  matrix.  However,  the  eigenvectors  of 
the  Zemike  covariance  matrix  serve  as  a  Karhunen-Loeve  basis  set.  These 
eigenvectors  can  represent  turbulence  because  they  are  not  correlated,  i.e. 
each  eigenvector  is  formed  by  superposition  of  Zemike  polynomials  in  such 
a  way  that  the  K,  are  orthogonal,  satisfying  the  first  objective.  Additionally,  the 
cofTesponding  eigenvalue  \  multiplied  by  (D/rj*^  (where  D  is  the  aperture 
diameter  and  r«  is  Fried's  cohererKe  lerHJth  (Fried,  1965))  gives  the  appropriate 
variarKe  for  tiuit  K,  spatial  compor^nt  corresponding  with  Kolmogorov 
turbulence,  satisfying  the  second  objective. 

Spedficalty,  the  low  spatial  frequerK:y  contribution  to  the  phase 
screen  can  be  expanded  in  terms  of  these  Karhuner)>Loeve  components 
(Cochrane,  1985) 

♦W  •  E  T, 

P»1 

where  p,^  represwits  the  rnimber  of  low  spatial  frec^jerKy  KartHJnen^.oeve 
terms  irKluded  and  the  coefticients  are  Gaussian  random  numbers  with 
variarKa  arn]  scaled  by  {Otrj^  to  the  specHtc  strer^gth  of  turbulecKe  used. 
The  simuiatkms  use  the  first  five  terms  ( p^  «  5  )  as  a  compromise  between 
completeness  of  tow  spatial  frequerKy  correction  and  oomputatkxMl  efflderKry. 
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Figure  (47)  shows  a  surface  plot  of  one  realization  of  the  first  two  terms  of  the 
correction,  which  are  very  close  to  x-  and  y-tilt  (i.e.  phase  terms  linear  in  x  and 
y,  respectively).  Figure  (48)  shows  a  realization  of  the  3rd,  4th,  and  5th 
correction  terms,  which  resemble  the  wavefront  aberrations  associated  with 
defocussing  and  astigmatism  in  a  conventional  imaging  system.  To  implement 
the  low  spatial  frequency  phase  correction,  (1)  the  scalar  product  was  formed 
between  the  initial  Founer  transform  phase  screen  and  each  Karhunen-Loeve 
function  K^,  jiving  the  relative  strength  of  that  in  the  initial  phase  screen;  (2) 
this  amount  of  each  spatial  component  was  then  subtracted  from  the  phase 
screen;  and  (3)  the  component  was  then  added  back  to  the  phase  screen  in 
the  proper  amount  given  by  the  product  of  a  Gaussian  random  number  with 
variance  \  and  the  factor  to  scale  to  ttie  particular  strength  of 
turbulence  used. 

Cochrane's  computer  routines  only  calculate  the  Karhunen-Loeve 
correction  terms  over  the  largest  ctrde  that  fits  inside  the  calculation  gnd,  as 
shown  in  Figs.  (47)  and  (48).  Correction  of  the  E-field  over  the  entire  computer 
simulation  grid  requires  the  Karhurien-Loeve  terms  be  calculated  over  an  area 
that  is  at  least  /?  larger  on  each  side  than  the  field  grid  (Cochrane,  1985).  If 
the  E-field  grid  size  N  is  chosen  as  a  power  of  2.  then  the  Karhunen-Loeve  grid 
must  be  2Nx2N.  These  simulations  however  used  only  an  NxN  grid 
(1024x1024)  because  the  Sun  SparcStations  could  not  handle  the  memory 
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Figm  47  Terms  1  and  2  of  Karttunen-Loeve  low  spatial  frequency  correction 
to  phase  screen,  represented  in  optical  path  length  for  wavelength  »  500  nm. 
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requirements  of  the  doubled  grid  (2048x2048)  without  a  major  revision  of  the 
propagation  code.  Since  the  source  hjnction  was  chosen  to  minimize  energy 
scattered  off  the  grid,  very  little  energy  fell  in  the  uncorrected  comers  of  the  grid 
and  the  simulation  remained  valid. 

Noll's  covariance  matrix  assumes  Kolmogorov  turbulence 
spectrum  with  zero  inner  scale  making  the  Karhunen-Loeve  correction  terms 
strictly  apply  only  to  this  spectrum.  These  simulations  however  incorporate 
nonzero  inner  scales  into  the  spectrum  of  refractive  index  fluctuations. 

Because  the  first  five  Karhunen-Loeve  terms  cover  scale  sizes  on  the  order  of 
the  grid  widtft,  which  is  much  larger  than  the  inner  scale  size  the  Karhunen- 
Loeve  correction  terms  derived  for  zero  inner  scale  remain  valid  for  correcting 
noruero  inner  scale  phase  screens. 

Cochrane  (1985)  showed  that  such  low  spatial  frequency 
correction  greatly  improved  the  phase  structure  function.  Two  terms  corrected 
the  structure  function  to  within  10%  of  the  theoretical  Kolmogorov  behavior,  and 
five  terms  corrected  to  within  5%.  compared  with  30  - 1000%  discrepancies 
without  any  correction  at  low  spatial  frequencies.  Figure  (49)  plots  the 
normalized  irradiance  variances  from  computer  simulations  as  a  function  of 
for  zero  Karhunen-Loeve  low  frequency  correction,  two-term  correction,  and 
five-term  correction.  Zero  correction  underestimated  the  irradiance  variance  in 
8)e  Rytov  regime  by  approximately  5%,  and  the  two-term  tilt-type  correction  also 


118 


Normolized  irradkince  vorione*  /  b«to  tquorsd 


Figura  49  Normalized  irradiance  variance,  normalized  by  versus 
turbulence  strength  for  zero  (dotted),  two^erm  (dashed),  and  five-term  (solid) 
Karhunen-Loeve  low  spatial  frequency  corrections  to  phase  screens. 


underestimated  by  about  5%.  The  five-term  correction  with  its  focussing  type 
terms  raised  the  variance  to  within  about  2%  of  the  theoretical  variance.  In  the 
saturation  regime,  the  two-  and  five-term  corrections  fit  better.  Figure  (50)  plots 
the  coherence  length  rg,  normalized  by  the  theoretical  coherence  length,  as  a 
function  of  turbulence  strength  for  zero,  two-,  and  five-term  Karhunen-Loeve  low 
spatial  frequency  corrections.  In  the  Rytov  regime  where  simulation  should 
dosely  approximate  theory,  the  zero  correction  overestimated  the  coherence 
length  by  approximately  35%.  The  tilt-type  correction  (two  terms)  estimated  the 
coherence  length  within  about  5%,  and  the  five-term  correction  achieved 
agreement  within  about  2%.  These  behaviors  formed  the  guideline  that  some 
type  of  low  spatial  frequency  correction  was  required  to  achieve  valid 
coherence  lengths  from  computer  simulation. 

This  low  spatial  frequency  correction  method  remained 
computationally  feasible  because  the  Noll  covariance  matrix  only  needed  to  be 
calculated  once  and  because  only  a  few  terms  of  the  Karhunen-Loeve 
expansion  were  used.  However,  the  code  can  become  memory  intensive 
because  it  saves  a  full  NxN  grid  for  each  Karhunen-Loeve  function  in  addition 
to  the  NxN  phase  screen  itself.  Implementation  with  2048x2048  or  larger  grids 
becomes  problematical  except  on  very  large  computers  with  about  one  gigabyte 
of  RAM.  Fried(1993)  has  proposed  a  simpler  x-  and  y-tilt  correction  algorithm 
and  indicates  that  this  performs  almost  as  well  as  including  the  higher  order 
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dotted  line  »  0  inner  scale,  0— term  KL  correction 
dashed  line  0  inner  seole.  2 -term  KL  correction 
solid  lino  ■  0  inner  scole,  5-term  KL  correction 


Flgiim  60  Spherical  wave  coherence  length,  normalized  by  theoretical 
coherence  length,  versus  turbulence  strength  for  zero  (dotted),  two-term 
(dashed),  and  five-term  (solid)  Karhunen-Loeve  tow  spatial  frequency 
corrections  to  phase  screens. 
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Karhunen-Loeve  terms,  while  being  less  memory  intensive  and  faster  than  the 
Karhunen-Loeve  correction  method.  Figs.  (49)  and  (50)  indicate  that  the  tilt- 
only  correction  may  systematically  underestimate  the  normalized  irradiance 
variance  by  about  5%  and  overestimate  the  coherence  length  of  a  spherically 
diverging  wave  by  about  5%  in  the  Rytov  regime.  Depending  on  the 
application,  the  Karhunen-Loeve  correction  to  higher  order  terms  may  prove  a 
useful  refinement. 


H.  COHERENCE  LENGTH 

The  coherence  length  of  the  E-field  was  calculated  from  the  atmospheric 
MTF  with  the  theory  outlined  in  Chapter  II.  but  the  actual  implementation  of  the 
MTF  calculation  and  parameterization  of  the  coherence  length  required 
careful  consideration.  The  calculation  methods  that  worlted  best  for  these 
investigations  are  discussed  in  this  section. 

Equation  (60)  of  Chapter  II  relates  the  atmospheric  MTF  and  the 
spherical  wave  structure  function  to  the  coherence  properties  of  the  E-field 


MTF^ 


m) 

9  * 


( W{r')  ) 


(134) 


Again,  U(f)  represents  the  E-field  and  W(f)  represents  the  aperture  function. 
The  autocorrelations  of  U  and  W  indicated  in  Eq.  (134)  could  be  done  by 
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summing  the  complex  products  of  the  E-fields  over  multiple  pairs  of  points  or  by 
implementing  the  autocorrelation  via  FFT  techniques.  Both  versions  were  tried 
for  these  investigations  and  produced  similar  results,  but  the  FFT  version 
provided  a  much  more  thorough  autocorretation  with  greater  computational 
efficiency. 

The  FFT  autocorrelation  technique  used  the  MCF  introduced  earlier. 


Equation  (57)  defined  the  mutual  coherence  function  as 

MCF  .  (  4/*(A.f,)  >,  (13S) 

which,  for  a  single  time  =  tj  =  t  and  assuming  homogeneity,  may  be  written 
as  the  spatial  autocorrelation  of  the  E-fieid 

MCF(F')  =  /  a  U(f  +  /")• 

Substituting  Eq.  (136)  into  the  Fourier  transform  identity, 

MCF{r')  =  f  I  dt"  MCF{r")  J  (137) 

gives 

MCF{r)  ^  I  dt[  I  dt  U{M")  I  e  I  (138) 

Rearranging  ffie  integrations, 

MCF{F)  ^  f  df[f  dt  U^if)  I  f  dt"  U{M")  e  )  J  (139) 

Changing  variables  to  ^  =  7  ♦  7",  7"  »  8  -  7,  dT "  » 
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MCF(f')  -  f  df  if  dt  t/*(0  I  /  di  U{i)  1  (140) 

The  inner  integral  is  the  Fourier  transform  of  U(l)  and  is  denoted  by  FTI  U  ]. 

The  next  innermost  integral  is  the  inverse  Fourier  transform  of  U*(/)  and  is 
denoted  by  IFT[  U*  1.  Then, 

MCF{F)  =  f  df  FTI  U]  tFT\U*]  (141) 

which  is  yet  another  inverse  Fourier  transform,  symbolically  written 

MCF{F)  -  IFT[  FT\U\  IFT\U*\  ).  (‘*^2) 

Equation  (142)  expresses  the  autocorrelation  (MCF)  of  the  E-field  in  terms  of 
Fourier  transform  techniques  easily  implemented  with  discrete  Fourier 
transforms. 

This  Fourier  transform  method  of  calculating  the  autocorrelation  of  a 
function  is  faster  and  more  complete  ttuin  the  more  laborious  technique  of 
averaging  the  products  of  the  E-field  at  point  pairs.  The  product  E*(/|  )E(f2  )  is 

complex,  but  because  of  the  averaging  that  occurs  in  the  autocorrelation,  tfre 
real  part  attains  a  stable  value  while  the  imaginary  part  averages  toward  zero. 

For  the  FFT  implementation  on  a  1024x1024  grid,  the  imaginary  part  is  typically 
*10**  while  the  real  part  ranges  between  0  and  numbers  on  the  order  of  unity. 
The  point  pairs  technique  produces  an  equivalent  real  part  but  reduces  the 
imagiriary  part  down  to  only  -'10'’.  These  investigations  used  the  FFT  version  for 
all  autocorrelations.  Additionally,  the  autocorrelations  were  normalized  by  their 
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zero  lag  values  to  ensure  that  the  magnitude  of  the  MTF,^  calculated  by 
Eq.  (134)  lies  between  0  and  1. 

Because  the  spherical  wave  structure  function  refers  to  points  across  a 
spherical  wavefront,  U  and  W  must  similarly  represent  the  E-field  across  a 
spherical  surface.  However,  the  computer  simulation  implements  the  E-field  on 
a  plane  perpendicular  to  the  axis  of  propagation  and  incorporates  die  spherical 
wave  nature  of  the  E-field  by  applying  a  quadratic  phase  curvature  across  the 
plane.  To  convert  from  this  plane  representation  of  the  E-field  in  the  simulation 
to  a  spherical  surface  representation  appropriate  for  the  autocorrelation 
calculations  of  Eq.  (134),  the  quadratic  phase  factor  across  the  plane  must  be 
removed  from  the  E-field.  This  effectively  assumes  that  the  amplitude  and  the 
phase  fluctuations  of  the  E-field  across  the  plane  of  the  computational  grid 
closely  approximate  the  E-field  across  the  spherical  wave.  This  assumption  is 
justified  because  the  maximum  physical  separation  of  the  true  spherical 
reference  surface  from  the  plane  surface  of  the  grid  is  at  most  1/50  the 
coherence  length  of  the  E-field.  However,  the  removal  of  the  quadratic  phase 
curvature  from  the  E-field  of  the  grid  pro\ms  crucial  to  the  autocorrelation  of  U 
because  the  quadratic  phase  ftictor  undergoes  ~130  multiples  of  2ic  phase 
change  between  the  center  and  the  outer  edge  of  the  grid  for  these  simulations 
and  would  otherwise  completely  obscure  the  actual  E-field  fluctuations. 
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Propagation  with  a  plane  or  beam  wave  requires  a  different  approach  to 
remove  the  proper  amount  of  phase  curvature.  For  a  plane  wave,  the 
wavefront  coincides  with  the  plane  of  the  grid  so  that  no  phase  curvature  needs 
to  be  removed.  A  pure  beam  wave  (Qaussian  profile)  exhibits  a  spherical 
phase  with  a  radius  of  curvature  larger  than  the  propagation  distance  L,  and 
this  phase  could  be  calculated  analytically  and  removed.  However,  the  use  of  a 
finite  E-field  confined  to  the  propagation  grid  introduces  phase  effects  other 
than  simple  quadratic  curvature.  Fortunately,  the  E-fieid  propagated  through 
zero  turbulence  contains  this  phase  curvature  information  (Walters,  1994).  For 
spherical  wave  propagations,  the  removal  of  the  phase  curvature  by  the 
analytical  calculation  and  by  using  the  zero  turbulence  propagated  field 
provided  identical  coherence  lengths.  These  investigations  used  the  latter 
method  for  the  coherence  length  calculations  for  all  beam-like  and  plane-wave- 
like  propagations. 

Equation  (134)  requires  that  the  autocorrelation  of  U  must  be  averaged 
over  multiple  realizations  to  achieve  the  long  exposure  MTF.  To  implement  this 
requirement,  the  autocorrelations  from  several  realizations  were  calculated  and 
averaged  together,  and  a  coherence  length  then  determined  via  Eq.  (134).  This 
method  produced  coherence  lengffis  that  agreed  with  theory  within  ->5%  in  the 
Rytov  regime.  Recalling  Fig.  (45)  for  ffie  coherence  lengff)  versus  number  of 
realizations  used  in  the  average,  the  number  of  fields  included  in  the  MTF 
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average  may  be  as  few  as  5,  though  20  provided  a  more  statistically 
reproducible  coherence  length.  Again,  these  simulations  used  20  realizations 
because  these  fields  had  already  been  generated  for  the  normalized  irradiance 
variance  calculations. 

Figure  (51)  plots  the  coherence  lengths  calculated  from  the  average  MTF 
and  also  plots  the  average  of  coherence  lengths  calculated  with  single 
realization  MTF's.  The  average  of  individual  realization  coherence  lengths 
exceeded  the  averaged  MTF  coherence  length  by  '>20%  at  low  turbulence 
strengths  but  eventually  agreed  within  5%  near  the  saturation  regime.  The 
single  realization  MTF's  showed  a  larger  coherence  length  at  low  turbulence 
strengths  because  the  contributions  frcm  low  spatial  frequency  components  had 
not  been  reduced  through  averaging.  At  these  low  turbulence  strengths, 
multiple  realizations  were  required  to  average  these  low  spatial  frequency 
contributions  and  to  achieve  the  appropriate  long  exposure  MTF.  At  higher 
turbulence  strengths  near  saturation,  the  E-field  had  more  energy  at  high  spatial 
frequencies  that  dominated  the  MTF.  The  autocorrelation  for  a  single 
realization  now  averaged  over  many  coherence  lengths  and  yielded  an  MTF 
dose  to  the  average  MTF  for  multiple  realizations. 

Using  a  finite  beam  to  approximate  a  spherically  diverging  wave 
produced  a  radial  dependence  of  the  average  irradiance  that  affected  the 
coherence  length  calculation,  but  only  to  a  small  degree.  To  compensate  for 
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mis  radial  deparKlenca,  each  E-field  was  divided  by  the  radial  average  E-field 
magnitude  from  the  30  realizations.  Similar  to  the  irradiance  variance 
calculation,  this  average  E-fieid  was  calculated  by  dividing  the  grid  into  rings 
one  grid  element  wide,  taking  the  square  root  of  the  average  intensity  for  each 
ling  over  the  30  realizations,  and  then  performing  an  area  weighted,  running 
mean  across  the  rings  to  smooth  the  variations.  While  this  radial  compensation 
proved  essential  for  the  normalized  irradiance  variance  calculation,  it  only 
changed  the  coherence  iengm  by  ~1%,  which  was  significantly  less  thnn  the  ~5% 
discrepancy  from  the  low  spatial  frequer>cy  correction. 

The  above  considerations  allowed  calculation  of  the  right-hand  side  of 
Eq.  (134)  for  the  atmospheric  MTF.  Because  of  the  statistical  nature  of  the 
propagation  through  turbulence,  no  set  of  realizations  yielded  an  atmospheric 
MTF  that  exactly  followed  the  exponential  rolloff  with  distance  predicted  by  the 
left-hand  side  of  Eq.  (134).  Methods  had  to  be  developed  to  parameterize 
these  atmospheric  MTF's  from  the  simulations  arid  extract  an  appropriate 
coherence  length  oorrespondirig  to  the  structure  tijnction  D(r). 

For  the  case  of  Kolmogorov  turbulerice,  Fried  (1966,  p.  1380-1383) 
derived  the  wave  structure  function  arid  expressed  it  in  terms  of  the  single 
cohererice  length  parameter  r, 

Oir)  >  6J8  (143) 
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where  n  =  5/3.  The  atmospheric  MTF  then  becomes 


iUTF, 


(144) 


To  extract  the  coherence  length  ro  arxt  the  exponent  n  (allowing  the  possibility 
that  n  may  vary),  take  the  natural  logarithm  of  both  sides  and  rearrange 
(Walters.  1993) 

,14.) 

Taking  the  natural  logarithm  again. 

ln(  -  ^  ) )  =  n  h(f)  -  n  ln(r,).  (14«» 

This  has  the  linear  form  y  «  ax  b,  where  y  represents  the  left  hand  side, 
a  =  n,  X  >  in(r),  and  b  =  -n  In(ro).  The  atmospheric  MTF  can  now  be 
characterized  with  two  parameters,  r,  and  n,  or  just  the  single  parameter  rg 
assuming  n  ■  5/3. 

To  im|[Nement  these  parameterizations.  the  atmospheric  MTF  was  first 
calculated  from  the  E-field  using  the  autocorrelation  methods  above  and  then 
radially  averaged  to  yield  MTF,g.^(r)  for  use  in  Eq.  (146).  A  linear  least 
squares  fit  calculation  provided  the  slope  a  -  n  and  the  intercept  n  In(ro),  giving 
rg.  To  obtain  the  single  parameter  characterization,  n  was  set  to  5/3  in 
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Eq.  (146)  and  least  squares  techniques  applied  to  obtain  the  intercept  and  thus 
the  single  parameter  Tq. 

Figure  (52)  shows  the  coherence  lengths  calculated  with  these  least 
squares  methods,  divided  by  the  Rytov  theory  coherence  length  for  Kolmogorov 
turbulence.  The  two  parameter  characterization  with  both  exponent  n  and 
coherence  length  rg  never  provided  a  consistent,  smoothly  varying  coherence 
length.  At  low  turbulence  strengths  in  the  Rytov  regime  where  the  E-fieid 
coherence  length  is  larger  than  the  calculation  aperture,  the  two  parameter  fit 
predicted  coherence  lengths  up  to  50%  higher  than  the  theoretical  value  and 
thus  appears  unreliable.  Additionally,  it  showed  an  anomalous  bump  around 
-  0.5  .  The  single  parameter  least  squares  technique  with  n  -  5/3 
accurately  characterized  the  coherence  length  within  5%  at  low  turbulence 
strength,  but  still  showed  the  bump  at  ~  0.5  . 

An  alternate  technique  to  obtain  a  single  parameter  characterization 
assumed  n  =  5/3  in  Eq.  (144)  and  used  a  binary-type  search,  or  iterative  fit,  to 
find  the  coherence  length  Tg  that  minimized  the  variance  between  the  average 
MTF,g^(r)  and  the  right-hand  side  of  Eq.  (144).  Though  not  analytical,  the 
resulting  Tg  gave  an  MTF  that  often  fit  the  actual  MTF,g^(r)  more  closely  by  eye 
than  the  least  squares  methods,  especially  for  low  turbulence  where  coherence 
lengths  were  larger  than  the  grid  size.  To  implement  the  technique,  an  initial 
large  range  of  Tg  (for  example,  0  to  100  m)  was  divided  in  half,  the  midpoint  rg's 
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Figura  S2  Coherence  lengths  calculated  by  one-parameter  (dashed  line)  and 
two-parameter  (dotted  line)  least  squares  techniques,  the  iterative  fit  technique 
(solid  line),  and  the  HWHM  (dash-dot  line). 
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of  each  half  weie  substituted  Into  Eq.  (144).  and  the  variances  were  calculated. 
Whichever  rg  provided  the  least  variance  became  the  middle  of  the  next  range, 
and  the  other  rg  became  the  new  high  or  low  boundary.  The  process  was  then 
repeated,  so  that  each  iteration  reduced  the  range  of  possible  rg  by  1A3.  This 
procedure  was  iterated  100  times  or  until  the  variance  was  less  than  1x10*'. 

Figure  (52)  shows  the  coherence  lengths  predicted  by  this  iterative  fit 
method  along  with  the  coherence  lengths  from  the  one*  and  two-parameter 
least  squares  fits,  all  divided  by  the  Rytov  theory  coherence  length  for 
Kolmogorov  turbulence.  The  iterative  fit  coherence  length  was  ^-5%  too  large  at 
low  turbulence  but  did  not  show  the  anomalous  bump  around  Pg^  ~  0.5  . 

In  the  saturation  regime  where  turbulence  is  high  and/or  path  lengths  are 
long,  multiple  scattering  becomes  significant  and  saturates  the  irradiance 
variance  (Martin  and  Flatty,  1988).  Correspondingly,  this  physical  phenomenon 
may  also  affect  the  coherence  length  and  modify  the  value  of  n  or  the  form  of 
Eq.  (144)  for  the  saturation  regime.  The  half  width  at  half  maximum  (HWHM) 
provided  a  coarser,  one  parameter  method  of  characterizing  the  atmospheric 
MTF  that  did  not  depend  on  any  assumptions  about  the  form  of  the  structure 
function  and  could  be  calculated  directly  from  the  atmospheric  MTF  by  linearly 
interpolating  between  the  pair  of  points  bounding  MTF,^„„,(r)  =  1/2.  Figure  (52) 
plots  the  HWHM,  divided  by  the  HWHM  predicted  by  Rytov  theory  assuming 
Kolmogorov  turbulence.  The  HWHM  plot  starts  at  Pg^  =  0.05  because  lower 
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turbulence  strengths  gave  a  coherence  targe  enough  that  the  MTF,^(r)  had 
not  reached  half  its  maximum  value  within  the  calculation  region.  The  HWHM 
lengths  followed  the  theoretical  values  ^thin  5%  in  the  Rytov  regime  and  did 
not  show  the  anomalous  bump  around  ~  0.5  .  Note  that  the  iterative  fit 
coherence  lengths  agreed  with  the  HWHM  plot  better  than  the  least  squares 
techniques.  The  HWHM  and  iterative  fit  proved  to  be  the  most  stable 
parameterizations  of  the  E-field  coherence  length,  and  were  used  in  all 
subsequent  coherence  length  comparisons  and  plots. 

Other  factors  were  also  considered  in  the  coherence  length  calculation. 
As  stated  earlier,  some  choices  for  point  source,  such  as  the  Airy-type  source 
required  more  energy  at  high  spatial  frequencies  than  others,  such  as  a  narrow 
Gaussian.  While  this  difference  produced  <  2%  effect  in  the  normalized 
irradiance  variance  calculations.  It  appeared  to  have  more  effect  on  coherence 
length  calculations.  Figure  (53)  shows  that  an  Airy-type  source  produced 
coherence  lengths  up  to  10%  higher  at  low  turbulence  strengths  than  those 
from  the  Gaussian-type  source  of  Eq.  (109).  Therefore,  only  the  latter  type 
source  was  used  for  further  investigations. 

For  an  arbitrary  spectrum  of  refractive  index  fluctuations  <1>„(k,z),  the 
integral  formulation  Eq.  (61)  provided  tiie  wave  structure  function  for  the 
atmospheric  MTF  in  Eq.  (134).  Spedfically,  numerical  integration  of  Eq.  (61)  for 
the  spectra  with  grid  cutoff,  Gaussian,  and  Hill/Frehlich  viscous  convective 
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Figura  53  Iterative  fit  coherence  lengths  versus  turbulence  strength  for  Airy- 
type  source  (dotted),  Gaussian-type  source  (solid),  and  Gaussian-type  source 
with  zero  Karhunen-Loeve  low  spatial  frequency  corrections  (dashed). 


enhancement  inner  scales  allowed  comparison  of  theory  (Eq.  (134))  and  the 
simulation  (described  in  Chapter  IV).  However,  in  the  structure 

function  numerical  integration,  care  had  to  be  exercised  in  properly  treating  the 
low  spatial  frequency  portion  of  the  integrand  since  the  spectra  contained  a 
Singularity  as  k  approached  zero.  This  low  frequency  portion  was  critical 
to  obtain  coherence  lengths  that  approached  the  Kolmogorov  theoretical  values 
in  the  Rytov  regime.  The  integral  was  successfully  evaluated  by  integrating 
analytically  for  0  <  k  <  ic^  .  ( where  1x10*^  m*' )  (Walters,  1994)  since  the 
inner  scale  function  F(k^)  -  1  here  and  this  portion  of  the  spectrum  remains 
Kolmogorov.  The  remaining  portion  of  the  integral  that  contained  the  inner 
scale  contribution  was  carried  out  numericaily.  More  realistic  spectra  could  also 
have  included  an  outer  scale,  but  again  no  universal  form  of  outer  scale  exists 
due  to  anisotropy  of  the  atmosphere  at  large  scale  sizes.  When  included  in  the 
numerical  integrations  for  test  purposes,  an  outer  scale  raised  the  coherence 
length  compared  to  the  Kolmogorov  case.  However,  these  investigations  did 
not  use  an  explicit  outer  scale  in  the  simulations. 
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IV.  RESULTS 


A.  NORMAUZED  IRRADIANCE  VARIANCE 

The  computer  simulation  guidelines  and  considerations  discussed  in 
Chaplir  III  were  implemented  to  investigate  the  behavior  of  the  normalized 
irradiance  variance  and  E-field  coherence  length  in  the  Rytov  and  saturation 
regimes  for  the  grid  cutoff,  Gaussian,  and  Hill/Frehlich  viscous-convective 
enhancement  inner  scales.  Specifically,  the  simulations  S4>ply  to  stratosphenc 
propagation  with  propagation  distance  L  =  200  km,  wavelength  X  -  500  nm. 
strengths  of  turbulence  =  [5x10**.  50],  and  inner  scale  sizes  [0,  15]  cm. 

The  simulations  used  a  1024x1024  grid  with  grid  element  size  given  by  Eq. 

(92),  an  Airy-type  source  modified  to  produce  a  final  zero  turbulence  irradiance 
pattern  with  edges  apodized  by  a  Gaussian  (Eq.  (109))  and  width  corresponding 
to  half  the  grid  width,  32  phase  screens  utilizing  a  five-term  Karhunen-Loeve 
low  spatial  frequency  correction,  and  30  realizations  in  each  set  of  runs.  The 
central  256x256  portion  of  each  propagated  E-field  was  used  for  the  normalized 
irradiance  variance  and  coherence  length  calculations. 

For  the  Gaussian  inner  scale  values  of  0  (grid  cutoff),  5, 10,  and  15  cm. 
Fig.  (54)  plots  the  normalized  irradiance  variance  versus  turbulence  strength 
in  the  Rytov  regime  from  both  numedcai  integration  of  the  equation  from  Rytov- 
Tatarski  theory,  Eq.  (34)  (dotted  lines),  and  from  computer  simulations 
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Figim  M  Normalized  irradiance  variance  (divided  by  from  simulation 
(solid)  and  numerical  integration  (dotted)  for  Gaussian  inner  scales  of  0,  5,  10, 
15  cm. 
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(solid  tines).  All  values  are  normalized  by  .  Numerical  integration  values 
for  0  cm  inner  scale  agreed  within  1%  of  the  theoretical  zero  inner  scale  values. 
The  difference  resulted  because  the  numerical  integration  was  limited  to  spatial 
frequencies  below  the  grid  cutoff  k^.  Larger  grid  cutoff  values  gave  closer 
agreement.  The  simulation  normalized  itradiance  variances  agreed  within  2% 
of  the  numerical  integration  values  for  all  four  inner  scale  values  examined. 
Nonzero  Gaussian  inner  scales  reduced  the  normalized  irradiance  variance 
below  the  zero  inner  scale  value  (by  10%.  25%,  and  40%  for  the  5,  10,  15  cm 
cases,  respectively).  Intuitively,  the  finite  inner  scale  suppressed  the  higher 
spatial  frequency  index  of  refraction  fluctuations  and  thus  reduced  the  variance. 

Figure  (55)  plots  the  normalized  irradiance  variance  (divided  by  Pg^)  for 
the  single  turbulence  strength  -  5x10**  and  Gaussian  Inner  scale  sizes  of  0 
(grid  cutoff),  5.  10,  and  15  cm.  The  numerical  integration  values  (dotted  line) 
and  computer  simulation  values  (solid  line)  showed  an  almost  linear  decrease 
of  the  normalized  irradiance  variance  with  increasing  inner  scale  size  in  the 
Rytov  regime.  As  Flatty,  Wang,  and  Martin  (1993)  point  out,  the  Gaussian 
inner  scale  does  not  accurately  describe  ttie  inner  scale  observed  in  the 
atmosphere  but  retains  usefulness  because  it  facilitates  some  theoretical 
calculations. 

For  the  Hill  viscous-convective  enhancement  inner  scale  sizes  of  0  (grid 
cutoff),  5, 10,  and  15  cm.  Fig.  (56)  plots  the  normalized  irradiance  variance  in 
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Flgun  65  Normalized  irradianoe  variance  from  simulation  (solid)  and  numerical 
integration  (dotted)  for  Gaussian  inner  scales  of  0,5,10.15  cm  at  s  5x10*^. 
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Figure  66  Normalized  irradiance  variance  (divided  by  from  simulation 
(solid)  and  numerical  integration  (dotted)  for  Hill  viscous-convective 
enhancement  inner  scales  of  0,5,10,15  cm. 
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the  Rytov  regime  versus  turbulence  strength  .  Numerical  integration  of  the 
Rytov-Tatarski  theory,  Eq.  (34)  (clotted  lines),  and  computer  simulation  (solid 
lines)  agreed  within  2%.  For  the  smaller  inner  scales,  the  normalized  irradiance 
variance  exceeded  the  zero  inner  scale  values  (by  30%  for  the  5  cm  case). 

The  viscous-convective  enhancement  of  tt)e  strength  of  higher  spatial 
vravenumber  fluctuations  near  the  inner  scale  wavenumber  increased  the 
variance.  Yet,  for  large  enough  inner  scale,  the  rolloff  beyond  the  enhancement 
suppressed  the  higher  spatial  frequency  fluctuations  enough  to  eventually 
reduce  the  variance  below  the  zero  inner  scale  variance  (by  30%  for  the  15  cm 
case).  The  10  cm  values  happened  to  lie  within  1%  of  the  0  cm  values. 

Figure  (57)  plots  the  normalized  irradiance  variance  (divided  by  for 
the  single  turbulence  strength  =  5x10**  and  the  Hill  viscous-convective 
enhancement  inner  scale  sizes  of  0  (gnd  cutoff),  2,  3.  4,  5,  6,  7,  10.  and  15  cm. 
Numerical  integration  of  the  Rytov-Tatarski  results  (dotted  line)  and  computer 
simulation  (solid  line)  dearly  illustrated  the  rising  and  then  falling  behavior  of 
the  normalized  irradiance  variance  wifri  increasing  inner  scale  size  in  the  Rytov 
regime.  The  normalized  irradiance  variance  achieved  a  maximum  for  ^  ~  4  cm, 
which  was  about  30%  of  the  Fresnel  length  R,  ■  l  /  i  n  -  12.6  cm. 

The  dashed  line  in  Fig.  (57)  diows  simulation  values  using  the  Frehlich 
parameterization  of  the  viscous-convective  enhancement  inner  scale.  The 
Frehlich  inner  scale  shifted  the  plot  slightly  to  smaller  inner  scale  sizes. 
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maximized  at  2%  less  than  the  Hill  maximum,  and  matched  the  Hill  values 
within  3%  over  the  range  of  inner  scale  plotted.  Additionally,  simulation  runs 
using  4  cm  Hill  and  Frehlich  inner  scales  agreed  within  3%  over  the  range  of 
turbulence  strengths  5x10**  <  <  50  .  Thus,  the  Hill  and  Frehlich  versions 

of  the  viscous-convective  enhancement  inner  scale  perform  almost  identically. 

Previous  investigations  have  iilusb'ated  the  dramatic  monotonic  rise  in 
normalized  irradiance  variance  in  the  saturation  regime  as  the  inner  scale  size 
increases  (Martin  and  Flatty,  1988).  Figure  (58)  shows  normalized  irradiance 
variance  from  computer  simulations  for  0  (grid  cutoff),  5,  10,  and  15  cm 
Gaussian  inner  scales  and  a  propagation  path  of  200  km.  Figure  (59)  shows  a 
corresponding  plot  for  0  (grid  cutoff),  5. 10,  and  15  cm  Hill  inner  scales.  The 
turbulence  values  range  from  the  Rytov  regime  (low  turbulence  with  1)  to 
the  saturation  regime  (high  turbulence,  sr  1  ).  The  Rytov  regime  showed 
again  the  behaviors  illustrated  with  Figs.  (54)  -  (57).  Increasing  Gaussian  inner 
scale  size  produced  monotonicaliy  decreasing  normalized  irradiance  variance. 
Martin  and  Flatt6  (1988,  1990)  provided  similar  plots  of  normalized  irradiance 
variance  with  a  Gaussian  inner  scale.  Tbe  Hill  viscous-convective 
enhancement  caused  the  normalized  irradiance  variance  to  rise  and  then  fall  as 
the  inner  scale  size  Increased.  However,  in  the  saturation  regime,  the 
normalized  irradiance  variance  increased  monotonicaliy  with  increasing  inner 
scale  size  for  both  the  Gaussian  and  Hilt  inner  scales.  The  transition  in 
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Flgura  58  Normalized  irradiance  variance  over  Rytov  and  saturation  regimes 
for  Gaussian  inner  scales  of  0  (solid),  5  (dashed),  10  (dash-dot),  and  15  cm 
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Figure  69  Normalized  irradiance  variance  over  Rytov  and  saturation  regimes 
for  Hill  viscous-convective  enhancement  inner  scales  of  0  (solid),  5  (dashed), 
10  (dash-dot),  and  15  cm  (dotted). 
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behavior  occurred  with  the  onset  of  saturation  around  ~  1  •  This  crossing 
behavior  has  been  plotted  for  the  log  intensity  variance  with  the  viscous- 
convective  inner  scale  by  Hill  and  Clifford  (1978). 

Figures  (54)  -  (57)  illustrate  the  close  agreement  between  numerical 
integration  and  computer  simulation  values  for  the  normalized  irradiance 
variance  at  low  strengths  of  turbulence.  This  agreement  provided  a  validity 
check  on  these  computer  simulations  that  incorporated  an  inner  scale. 

B.  COHERENCE  LENGTH 

Coherence  lengths  of  the  E-field  were  calculated  from  the  same 
1024x1024  simulation  runs  used  for  the  normalized  irradiance  variance 
calculation.  The  average  atmospheric  MTF  was  formed  from  30  realizations 
using  the  FFT  autocorrelation  method,  and  then  the  corresponding  coherence 
length  rg  and  HWHM  were  calculated.  The  iterative  fit  rg's  were  used  for 
comparisons  because  they  most  closely  followed  the  HWHM  behavior.  The 
HWHM  may  provide  the  coarsest  measure  of  coherence  length  but,  since  it 
requires  no  assumptions  about  the  form  of  the  MTF,  it  may  also  be  the  most 
reliable. 

Figure  (60)  shows  the  coherence  length  Tq  from  numerical  integration  of 
Eq.  (61)  for  an  approximately  zero  inner  scale  and  normalized  by  the  theoretical 
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Figure  60  Coherence  length  from  numerical  integration  for  zero  inner  scale 
versus  turbulence  strength.  Values  normalized  by  Kolmogorov  turbulence  zero 
inner  scale  coherence  length. 
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Kolmogorov  turbulence  zero  inner  scale  coherence  length  rg  of  Eq.  (66).  This 
plot  indicates  that  the  numerical  integration  coherence  lengttis  calculated  in 
these  investigations  were  accurate  within  2%. 

Figure  (61)  shows  the  coherence  length  rg  from  numerical  integration  of 
Eq.  (61)  with  the  grid  cutoff  =  318  rad/m  for  the  1024x1024  grid.  The  grid 
cutoff  did  not  affect  the  coherence  lengtti  until  «  1.5.  Above  that  level  of 
turbulence,  the  absence  of  the  very  high  spatial  frequency  contribution  to  the 
integral  caused  the  coherence  length  to  become  larger  than  the  theoretical  zero 
inner  scale  value.  The  grid  cutoff  inner  scale  implicit  in  the  computer 
simulations  with  the  finite  grid  should  have  a  similar  effect  at  these  high 
turbulence  values. 

Figure  (62)  shows  the  coherence  length  rg,  normalized  by  the  theoretical 
coherence  length,  from  computer  simulation  of  a  spherically  diverging  E-field. 

In  the  Rytov  regime,  the  simulation  agreed  with  theory  to  within  the  ~5% 
overestimation  due  to  using  only  five  terms  in  the  Karhunen-Loeve  low  spatial 
frequency  correction  to  the  phase  screens  (see  Fig.  (50)).  In  the  saturation 
rr  ime,  the  simulation  coherence  length  (solid  line)  dropped  below  the 
theoretical  prediction  (by  ~25%  at  Pg^  =  50)  and  the  HWHM  (dotted  line)  mirrored 
this  decrease.  For  strong  turbulence,  the  first  order  perturbation  theory  basis 
for  coherence  length  appears  to  lose  validity,  as  it  did  for  the  normalized 
irradiance  variance  in  the  saturation  regime. 


149 


0.0001  0.0010  0.0100  0.1000  1.0000  10.0000  100.0000 

Beta  squared 


FIgura  Coherence  length  from  numerical  integration  for  grid  cutoff  inner 
scale  versus  turbulence  strength.  Values  normalized  by  Kolmogorov  turbulence 
zero  inner  scale  coherence  length. 
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Flgura  62  Coherence  length  r,  (solid  line)  and  HWHM  (dotted  line)  for 
spherical  wave  propagation  through  turbulence  (Values  normalized  by 
theoretical  coherence  lengths). 


The  drop  in  coherence  length  presumably  came  from  strong  scattering  in 
high  turbulence,  but  a  possible  origin  in  the  computer  simulation  was  also 
considered.  The  simulations  started  witr.  an  approximate  point  source  to 
emulate  spherical  divergence  of  the  E-fieid  and  required  the  E-field  to  remain 
basically  confined  within  the  simulation  grid  over  the  propagation.  The  resulting 
E-field  properties  could  have  differed  from  those  of  a  true  spherical  wave.  To 
investigate  this  possibility,  the  width  of  tfie  final  irradiance  pattern  was  varied 
such  that  wider  final  irradiance  fields  (hence  narrower  sources)  more  dosely 
approximated  a  true  point  source.  Figure  (63)  plots  the  results  and  shows  that 
increasing  the  final  irradiance  width  (diamond,  then  triangle,  then  square,  then 
X)  actually  lowered  the  coherence  length  in  the  saturation  regime  while  still 
following  the  theory  in  the  Rytov  regime.  This  behavior  indicated  that  the 
observed  decrease  in  the  saturation  regime  was  physical  and  not  due  to 
simulation  constraints. 

To  further  investigate  this  phenomenon,  beam  wave  (i.e.  divergence 
intermediate  between  spherical  and  plane  waves)  and  plane  wave 
approximations  were  propagated  on  a  512x512  grid  in  which  the  width  of  the 
source  varied  from  4  grid  elements  (Ax)  to  364  grid  elements  (3/4  of  the  grid 
width).  Figure  (64)  plots  the  resulting  coherence  lengths  and  indicates  that, 
based  on  the  behavior  in  the  Rytov  regime,  the  4  Ax  source  (circles)  produced 
the  spherical  wave  having  large  divergence  of  the  E-field,  the  intermediate 
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FIgm  63  Coherence  lengths  for  varying  amounts  of  spherical  divergence, 
hence  final  irradiance  pattern  width:  diamond  4/8,  triangle  »  5/8,  square  =  6/8, 
X  =  7/8  grid  width. 
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Figim  64  Coherence  lengths  for  increasing  source  width  on  a  512x512  grid: 
cirdess4Ax  (spherical  wave);  X's^BAx,  squaressIGAx.  triangles=32Ax, 
diamonds=64Ax,  dotss128Ax.  asterisks=256^,  pluses=384Ax  (plane  wave). 
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source  sizes  (X^S  Ax,  square  »  16  Ax,  triangle  =  32  Ax,  diamond  »  64  Ax) 
produced  beam  wave  divergences,  and  wide  sources  (dot  -  126  Ax,  asterisk  = 
256  Ax,  plus  =  384  Ax)  produced  plane  waves. 

However,  as  the  turbulence  strength  approached  the  saturation  regime, 
the  behaviors  changed.  The  spherical-type  propagation  coherence  length 
dropped  ~15%  at  =  15,  the  beam  wave  propagations  actually  increased  in 
coherence  length,  and  the  plane  wave  propagations  first  decreased  ~5%  before 
increasing  ~15%  at  =  15.  Some  of  the  unevenness  in  the  beam  wave 
coherence  lengths  occurred  because  smaller  regions  of  the  E-field  were  used  to 
calculated  the  coherence  lengths  due  to  the  relatively  small  size  and  divergence 
of  these  waves. 

The  cause  of  these  behaviors  requires  hjrther  investigation,  but  a 
hypothesis  can  be  made.  As  noted  earlier,  the  spherical  wave  coherence 
length  probably  decreased  below  theory  for  strong  turbulence  where  the  Rytov- 
Tatarski  first  order  perturbation  theory  was  no  longer  valid.  Strong  scattering 
may  have  induced  spherical  divergence  of  the  beam  waves,  increasing  the 
coherence  length,  and  caused  the  plane  wave  approximations  to  diverge  like 
beam  waves  for  high  turbulence. 

Recapping  the  above  results,  investigations  of  coherence  length  via 
computer  simulation  indicated  that  first  order  perturbation  theory  for  coherence 
lengths  loses  validity  in  the  saturation  regime,  just  as  it  did  for  normalized 
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irradiance  variance.  The  behavior  of  coherence  length  in  the  saturation  regime 
for  spherical,  beam,  and  plane  wave  cases  requires  further  research. 

These  investigations  then  examined  the  effect  of  inner  scale  upon 
coherence  length  for  a  spherically  diverging  beam,  utilizing  the  1024x1024 
realizations  run  for  the  normalized  irradiance  variance  calculations.  Figure  (65) 
shows  the  coherence  length  r^  from  numerical  integration  of  Eq.  (61)  for 
Gaussian  inner  scales  of  5,  10,  and  15  cm.  The  Gaussian  inner  scale  made 
the  coherence  lengths  larger  by  reducing  the  energy  at  high  spatial  frequencies. 
Larger  inner  scales  monotonically  produced  larger  coherence  lengths  from 
numerical  integration  at  a  given  turbulence  strength  (~8%  larger  at  ~  1 
(beginning  of  saturation  regime)  for  ^  =  15  cm).  A  plot  of  HWHM  from 
numerical  integration  for  Gaussian  inner  scales  would  appear  similar  since  the 
theoretical  HWHM  is  proportional  to  the  Kolmogorov  coherence  length  r,, 

HWHM  =  0.382  r^.  (147) 

Figure  (66)  shows  the  coherence  length  r^  and  Fig.  (67)  shows  the 
HWHM  from  wave  optics  computer  simulations  with  Gaussian  inner  scales  of  5, 
10,  and  15  cm  (solid  lines),  superimposed  upon  the  numerical  integration 
predictions  of  Fig.  (65)  (dotted  lines).  ITie  rg  and  HWHM  simulation  values 
agreed  well  with  theory  for  <  0.1,  but  started  deviating  frc..i  theory  even 
before  Po^  =  1 ,  i.e.  the  onset  of  saturation  for  the  normalized  irradiance 
variance.  In  the  saturation  regime,  computer  simulation  r^  and  HWHM  still 
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Figm  66  Coherence  length  r,  from  numerical  integration  for  Gaussian  inner 
scales  of  5,  10,  and  15  cm,  normalized  by  the  theoretical  zero  inner  scale 
coherence  lengths. 
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Figura  66  Coherence  lengths  Tq  from  computer  simulation  with  Gaussian  inner 
scales  of  5,  10,  and  15  cm;  values  normalized  by  the  theoretical  zero  inner 
scale  coherence  lengths. 
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Figure  67  HWHM's  from  computer  simulation  with  Gaussian  inner  scales  of  5, 
10,  and  15  cm;  values  normalized  by  the  theoretical  zero  inner  scale  HWHM's. 
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agreed  well  with  each  other,  and  their  behavior  showed  a  combination  of 
decrease  in  coherence  length  due  to  saturation  regime  turbulence  and  inaease 
in  coherence  length  due  to  inner  scale. 

To  elucidate  the  effects  of  inner  scale  by  itself,  Figs.  (68)  and  (69)  show 
the  same  plots  of  computer  simulation  coherence  length  and  HWHM,  but  now 
normalized  by  the  computer  simulation  zero  inner  scale  and  HWHM  values, 
respectively,  effectively  removiny  .. j>aturation  contribution.  These  plots 
clearly  show  that,  even  in  the  saturation  regime,  inner  scale  increased  the 
coherence  length  (solid  lines)  similarly  to  the  predictions  of  theory  (dotted  lines). 

The  next  set  of  figures  plots  coherence  length  and  HWHM  for 
propagations  through  turbulence  with  the  Hill  viscous-convective  enhancement 
inner  scale.  Figure  (70)  plots  the  predicted  coherence  lengths  from  numerical 
integration  of  Eq.  (61).  Note  that  the  numerical  integration  coherence  lengths 
dropped  ~5%  in  the  range  =  [0.1,  1]  before  increasing  for  higher  turbulence 
strength.  Figures  (71)  and  (72)  plot  the  computer  simulation  coherence  lengths 
Tq  and  HWHM,  normalized  by  the  theoretical  values  for  spherical  waves.  In  the 
Rytov  regime,  the  coherence  length  r^  decreased  (~5%  for  ^  =  15  cm)  even  for 
very  low  turbulence  (Og^  <  0.1  )z,  but  then  increased  in  the  saturation  regime  as 
for  the  Gaussian  inner  scale.  Figures  (73)  and  (74)  plot  the  same  coherence 
lengths  r^  and  HWHM's  but  divided  by  the  zero  inner  scale  computer  simulation 
values  to  emphasize  the  effect  of  inner  scale.  Figure  (73)  clearly  shows  the 


small  decrease  in  coherence  length  r,  in  the  Rytov  regime,  and  both  plots  show 
the  general  agreement  between  theory  and  computer  simulation  for  the  effect  of 
inner  scale  upon  coherence  length  in  the  saturation  regime.  In  summary,  the 
inner  scale  increased  the  coherence  length  in  the  saturation  regime  as  much  as 
50%  compared  to  the  zero  inner  scale  case,  and  more  than  compensated  for 
the  decrease  from  strong  scatter. 
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Figure  68  Coherence  lengths  from  computer  simulation  for  Gaussian  inner 
scales  of  5,  10,  and  15  cm;  values  normalized  by  computer  simulation  zero 
inner  scale  rg. 
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Figure  69  HWHM's  from  computer  simulation  with  Gaussian  inner  scales  of  5, 
10,  and  15  cm;  valued  normalized  by  computer  simulation  zero  inner  scale 
HWHM’s. 
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Figim  70  Coherence  length  r,  from  numerical  integration  for  Hill  viscous- 
convective  enhancement  inner  scales  of  5, 10,  and  15  cm,  normalized  by  the 
theoretical  zero  inner  scale  coherence  length. 
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FIgtm  71  Coherence  lengths  r,  from  computer  simulation  with  Hill  viscous- 
convective  enhancement  inner  scales  of  5.  10,  and  15  cm;  values  normalized 
by  the  theoretical  zero  inner  scale  coherence  lengths. 
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Figura  72  HWHM's  from  computer  simulation  with  Hill  viscous-convective 
enhancement  inner  scales  of  5. 10,  and  15  cm;  values  normalized  by  the 
theoretical  zero  inner  scale  HWHM's. 
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Flgiire  73  Coherence  lengths  from  computer  simulation  for  Hill  viscous- 
convective  enhancement  inner  scales  of  5. 10,  and  15  cm;  values  normalized 
by  computer  simulation  zero  inner  scale  rg. 
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Figure  74  HWHM's  from  computer  simulation  with  Hill  viscous-convective 
enhancement  inner  scales  of  5,  10,  and  15  cm;  valued  normalized  by  computer 
simulation  zero  inner  scale  HWHM's. 
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V.  CONCLUSIONS 


Variations  in  the  index  of  refraction  within  a  turbulent  medium  alter  an 
E-field  propagating  through  the  medium.  The  Rytov-Tatarski,  perturbation 
theory  predicts  the  effects  of  the  turbulence  upon  the  irradiance  statistics  and 
coherence  length  of  the  propagating  E-fieid.  Computer  simulations  modeled  the 
propagation  of  the  E-field  through  the  turbulent  medium,  producing  irradiance 
and  coherence  statistics  to  compare  with  theoretical  results. 

These  investigations  used  a  split-step.  Huygens-Fresnei,  wave  optics, 
computer  simulation  to  model  an  E-field  propagating  through  a  turbulent 
stratosphere.  The  limits  of  validity  of  the  simulations  were  determined  based 
upon  aliasing  considerations,  choice  of  source,  and  robustness  of  statistical 
calculations,  and  produced  the  following  guidelines; 

•  The  element  size  for  an  NxN  grid  should  satisfy  Ax  =  l  /  K  . 

•  The  maximum  turbulence  strength  for  which  an  NxN  grid  produces  valid 
E-fields  is  given  by  *  0.1  N  . 

•  A  coherence  length  rg »  2.5  Ax  corresponds  to  the  maximum  turbulence 
strength  for  which  an  NxN  grid  produces  valid  E-fields. 

•  The  number  of  phase  screens  should  be  ^  30  . 

•  The  number  of  realizations  should  be  ^  30  . 

•  Low  spatial  frequency  corrections  to  phase  screens  improve  the  accuracy 
of  the  normalized  irradiance  variance  by  5%  and  of  the  E-field  coherence 
length  by  30%. 
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•  Half  width  at  half  maximum  of  the  atmospheric  MTF  and  an  iterative  fit 
provide  the  most  stable  parameterizations  of  the  E-field  coherence  length. 

•  Telltale  signs  of  aliasing  include  a  fine-grained  irradiance  pattern,  a  boxed 
perimeter  of  the  irradiance  pattern,  and  peaking  of  energy  toward  the 
center  of  the  computation  grid. 

This  research  investigated  the  effect  of  (1)  zero  inner  scale,  (2)  Gaussian 
inner  scale,  (3)  Hill's  and  (4)  Frehlich's  viscous-convective  enhancement  inner 
scales,  and  (5)  grid  cutoff  inner  scale  on  the  normalized  irradiance  variance  of  a 
spherical  wave  propagating  through  a  turbulent  medium.  For  the  Rytov  regime, 
the  normalized  irradiance  variances  grid  cutoff,  Gaussian,  and  Hill/Frehlich 
viscous-convective  enhancement  inner  scales  were  compared  to  the  zero  inner 
scale  case  and  to  the  values  from  numencal  integration  of  the  Rytov-Tatarski 
predictions.  For  low  turbulence  strengths,  the  variances  obtained  from  the 
simulations  agreed  within  2%  of  the  values  from  the  numerical  integrations. 

The  grid  cutoff  inner  scale,  implicit  in  discrete  grid  wave  optics  computer 
simulations,  affected  the  variance  negligibly  compared  to  a  true  zero  inner  scale 
at  low  turbulence  strengths  with  the  large  1024x1024  grid.  Application  of  a 
Gaussian  inner  scale  reduced  the  normalized  irradiance  variance  as  much  as 
40%  for  small  compared  to  the  zero  inner  scale  case.  The  more  realistic 
Hill  viscous-convective  enhancement  inner  scale  raised  the  normalized 
irradiance  variance  by  up  to  30%  for  smaller  inner  scale  values,  but  for  larger 
values  of  inner  scale  eventually  reduced  the  variance  below  the  zero  inner 
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scale  value  as  much  as  30%.  The  latter  contrasted  with  the  behavior  in  the 
saturation  regime  (f)o^  >  1)  where  larger  inner  scales  continually  enhanced  the 
normalized  irradiance  variance.  The  Frehii(^  parameterization  of  the  viscous- 
convective  enhancement  gave  normalized  irradiance  values  that  agreed  within 
3%  of  the  Hill  inner  scale  values  over  the  entire  range  of  turbulence  strengths 
investigated. 

The  coherence  of  the  E-field  was  studied  by  computing  the  average 
atmospheric  MTF  from  the  propagated  fields.  Parameterizing  the  MTF  with  a 
coherence  length  r,,  and  half  width  at  half  max  (HWHM)  allowed  comparison  of 
the  coherence  length  of  the  E-field  with  the  predicted  coherence  length  from  the 
Rytov-Tatarski-Fried  theory.  In  the  Rytov  regime,  simulation  coherence  lengths 
and  HWHM's  for  spherical  and  plane  wave  approximations  agreed  within  5%  of 
the  theoretical  coherence  lengths/HWHM’s  for  zero  inner  scale.  However,  in 
the  saturation  regime,  the  spherical  wave  coherence  length  decreased  as  much 
as  25%  below  the  theory.  Similar  decreases  resulted  for  different  widths  of  the 
final  E-field.  Beam  wave  approximations  gave  coherence  lengths  that 
increased  toward  the  spherical  wave  values  in  the  saturation  regime,  while 
plane  wave  approximations  deviated  from  ~5%  below  to  ~15%  above  the  theory 
in  the  saturation  regime.  Increasing  inner  scale  increased  the  coherence  length 
of  a  spherically  diverging  E-field  by  up  to  50%  relative  to  the  zero  inner  scale 
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case  in  the  saturation  regime.  The  amount  of  the  increase  agreed  with 
numerical  predictions  from  the  analytic  theory. 

Several  avenues  for  further  research  exist.  Larger  grid  sizes  and  finer 
mesh  could  explore  the  behavior  of  normalized  irradiance  variance  and 
coherence  length  at  higher  turbulence  strengths.  In  particular,  the  behavior  of 
coherence  length  for  spherical,  beam,  and  plane  waves  in  strong  turbulence 
conditions  need  a  more  precise  description.  The  atmospheric  MTF's  from 
simulation  and  theory  could  be  compared  directly,  rather  than  through  a 
coherenca  length  parameterization,  and  the  structure  functions  could  also  be 
calculated  and  compared  directly  to  investigate  whether  the  2/3  power  law 
structure  function  holds  in  saturation.  Simulations  could  allow  C„\z)  to  vary 
along  the  path  and  could  include  an  outer  scale  in  the  spectrum  ot  refractive 
index  fluctuations  to  study  the  effects  of  large  scale  anisotropy  of  atmospheric 
turbulence  on  normalized  irradiance  variance  and  coherence  length  (possibly 
basing  the  C„^(z)  variations  upon  high  frequency  radar  data  that  can  resolve  the 
large  scale  variations  down  to  about  300  m).  Larger  (  >  1  Gbyte  RAM)  and 
faster  computing  resources  would  provide  the  catalyst  for  all  such  further 
investigations  of  the  propagation  of  an  E-field  through  turbulence. 
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APPENDIX 


Sample  listing  of  YAPS  event  input  file: 


i>^NTS: 


BLPI^NAJIQN: 


t  7  ;clebug  flag  and  random  number  seed 

0.0  0.053033  izenith  angle  and  r_0  at  0.5  microns 

2  ;number  of  wavelengths 

0.5e-6  0.5e-6  ;list  of  wavelengths 

'surf  'beam'  0  0  0  2.531058298  0.0  0.0  0.00 


;surf  name,  vertex  location,  clear  aperture, 
;super-Gaussian  exponent,  inner  scale 

'surf  'atmosr  0  0  6250  8.0  0.0  0.0  0.00 


;surf  name,  vertex  location,  dear  aperture, 
;super-Gaussian  exponent,  inner  scale 
'prof  1024  1024  0.009886946  'none'  'dummy' 

;surface  size,  grid  element  size,  file  flags 
'end'  ;end  of  surface  summary 

2  1024  1024  ;number  of  fields  and  dimensions 

'times'  0.0  0.090002  ;time  initialization 

'thread'  0.0  1.0  1  .propagation  start 

finit'  1  0.009886946  2  0  0  0  1  0  0  200000  -i-l 


'apsrf  1  1  'beam' 

'aptou'  1 

'chgfcs'  1  0  0  1.0e+30 
'prop'  1  0  0  200000 
•fldcp'  1  2 
'prop'  2  0  0  0 

'openfl'  7data/davis/c118'  11 
'chgfcs'  2  0  0  200000 
'svflddx'  2  1 1  385  640  385  640 
fldcp'  1  2 


;initialize  field 

;apply  aperture  profile 

jconvert  from  amplitude/phase  to  complex 

jchange  focus  to  apply  spherical  phase 

;back  propagate  to  create  source 

;copy  field  to  use  as  source  later 

ipropagate 

;open  field  output  file 

;remove  spherical  phase  from  field 

;save  field  to  output  file 

;copy  source  field  to  working  grid 
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( following  three  steps  repeated  32  times  to  propagate  a  distance  L  ) 

'prop'  2  0  0  193750  ;propagate  field 

'mkscm'  6250.0  1.0e-18  193750.0  'atmosi' 

icreate  phase  screen  for  Az.  Cn2,  position 
'apsrf  2  1  'atmosi'  :apply  phase  screen  to  field 

(  save  the  field  ) 

'chgfcs'  2  0  0  200000  ;remove  spherical  phase  from  field 

'svfiddx'  2  11  385  640  385  640  ;save  field  to  output  file 

fidcp'  1  2  '.copy  source  field  to  working  grid 

(  repeat  above  propagation  30  times  } 


'closefi'  1 1  idose  output  file 

'end'  ;end  simulation 
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